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Implementation  Report 

In  general,  highway  and  airfield  pavements  consist  of  a  system  of  layers  resting 
on  a  foundation.  In  situ,  the  layers  can  be  nonlinear  and  vary  in  thickness,  homogeneity 
and  isotropy.  They  have  boundary  and  interaction  conditions  that  significantly  affect 
their  response  and  performance.  Loading  can  be  static,  moving  or  dynamic.  In  spite  of 
these  conditions  the  pavement  industry  has,  by  use,  largely  defined  rational  pavement 
analysis  and  design  as  a  linear  or  stepwise  linear,  layered  elastic  problem.  In  doing  so, 
the  layers  are  assumed  to  be  linear,  homogeneous  and  isotropic.  The  layers  are  assumed 
to  extend  to  infinity  and  interactions  are  limited.  Loading  is  taken  to  be  static. 

The  major  reason  given  for  continued  emphasis  on  use  of  the  layered  elastic 
model  is  that  it  is  easy  to  use.  There  are  various  layered  elastic  programs  that  run  on 
personal  computers  (PC).  Even  on  a  PC,  solutions  are  obtained  in  a  fraction  of  a  second. 
Models  based  on  the  finite  element  method  of  analysis  can  more  realistically  represent 
complex  pavement  problems.  Industry  has  resisted  use  of  such  models,  citing 
complicated  input  and  need  for  workstation  or  mainframe  computers.  Also,  computing 
times  are  large.  In  fact,  PCs  have  become  much  more  powerful  and  alternate  algorithms 
for  conducting  FEM  analyses  have  been  developed. 

A  study  was  conducted  that  involved  development  of  a  FEM  algorithm  that  would 
run  on  a  PC.  The  algorithm  is  an  explicit  solution  and  involves  a  vector  formulation 
representing  the  equations  of  motion.  Both  two-  and  three-dimensional  versions  of  the 
program  have  been  developed.  The  program  is  in  a  modular  format.  Initial  libraries  of 
properties  have  been  provided  for  the  load  and  material  modules,  respectively.  Additions 
can  be  made  to  the  libraries  as  needed.  Loads  can  range  from  static  to  dynamic.  Material 
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models  include  linear  and  nonlinear  elastic,  viscoelastic  and  plastic  or  a  combination. 
Both  the  loads  and  material  models  have  been  verified. 

The  three-dimensional  program  has  been  used  to  predict  deflection  response  of 
jointed  concrete  pavement  subjected  to  an  impulse  load  from  a  falling  weight 
deflectometer  (FWD).  The  predicted  response  was  in  close  agreement  with  in  situ 
measurements.  Other  example  problems  have  been  demonstrated  showing  that  the 
programs  can  be  utilized  for  the  general  analyses  of  pavement  systems  (concrete  and 
asphalt)  subjected  to  transient  dynamic  loading  histories. 

The  study  advisory  committee  recommended  that  the  report  and  programs  be 
accepted  and  therefore  satisfy  project  goals.  As  a  preliminary  of  the  committee  arriving 
at  this  recommendation,  a  presentation  was  made  on  options  for  developing  input  files 
and  displaying  results.  The  committee  decided  that  general  implementation  would  be 
greatly  enhanced  by  developing  a  graphical  interface  for  both  input  and  output. 
Consequently,  the  committee  recommended  that  an  implementation  study  proposal  be 
prepared  with  the  goal  of  developing  the  desired  graphical  interface. 


CHAPTER  1.  INTRODUCTION  AND  LITERATURE  SURVEY 

The  basis  of  a  rational  pavement  design  procedure  is  the  analysis  of  a  layered  pavement 
system.  This  pavement  analysis  system  should  include  models  of  the  solid  media,  material 
models  for  the  pavement  and  subgrades,  kinematics  conditions,  and  more  importantly, 
realistic  simulations  of  the  loading  conditions.  Considerable  efforts  have  been  reported  in 
recent  years  which  use  a  linear  elastic  analysis  to  obtain  the  basic  informations  needed  for 
the  development  of  such  a  design  procedure.  In  the  analysis,  the  pavements  are  assumed 
to  be  linearly  or  bi-linearly  elastic,  homogeneous,  and  isotropic.  Loading  are  usually  taken 
to  be  static.  As  these  highly  simplified  assumptions  do  not  correspond  to  the  real 
pavement  properties,  many  cases  exist  where  measured  pavement  responses  in  term  of 
deflection  or  strain  do  not  agree  with  those  predicted  by  a  linearly  elastic  model.   Shift 
factors,  for  example,  have  to  be  introduced  to  compensate  the  differences.  These  factors 
may  have  a  wide  range  in  magnitude  which  make  their  reliability  and  accuracy  susceptible 
to  large  errors. 

Continued  progress  toward  rational  pavement  design  that  would  effectively  simulate 
various  types  of  static  and  dynamic  loadings,  and  properties  of  conventional  and  emerging 
paving  materials  suggests  that  an  efficient  three  dimensional  finite  element  computer 
program  would  be  highly  desirable.  In  the  past,  such  programs  were  restrictive  and 
complicated,  and  required  the  use  of  large  mainframe  computer  systems.  The  advent  of 
new  algorithms  and  hardwares  in  recent  years  indicate  that  such  an  analysis  can  be 
performed  efficiently  on  a  personal  computer  system.  The  primary  motivation  of  this 
project  is  then  to  implement  such  an  efficient  computational  analysis  system,  which  may 
serve  as  the  foundation  of  a  rational  pavement  design  procedure. 


1 . 1  Explicit  Approach  of  Finite  Element  Analysis 

In  this  project,  an  explicit  approach  of  finite  element  analysis  is  adopted  as  the  basic 
methodology  for  the  development  of  a  pavement  analysis  system.  The  basic  ingredients  of 
this  approach  consist  of  a  vector  formulation  of  finite  elements  to  obtain  a  set  of  equations 
of  motion,  an  explicit  time-integration  procedure  to  find  solutions  of  the  equations  of 
motion,  and  a  co-rotational  formulation  to  handle  large  rotations. 

a.  A  Vector  Formulation  of  Finite  Elements: 

In  the  traditional  finite  element  analysis,  the  stiffness  and  force  matrices  are  calculated 
based  on  structural  discretization.  These  matrices  are  assembled  to  form  a  system  of 
simultaneous  linear  equations  for  the  solution  of  nodal  displacements.  For  dynamic 
problems,  a  similar  procedure  can  be  used  to  find  the  mass  matrix,  and  a  system  of 
differential  equations  of  motion  are  obtained. 

Instead  of  using  the  matrix  formulation,  a  transient  formulation  developed  earlier  by  Key, 
Belytschko  and  Hallquist  was  adopted.  In  the  formulation,  the  continuous  medium  is 
approximated  by  discrete  mass  particles.  Using  the  standard  finite  element  analysis,  energy 
equivalent  internal  and  external  forces  are  found.  They  are  the  forces  applied  on  the  mass 
particles.  Newton's  law  of  motion  and  time  integrations  complete  the  formulation  to 
determine  the  acceleration,  velocity  and  the  displacement  of  each  particle  for  a  particular 
time  increment. 

Since  the  discretization  considers  only  the  lumped  mass  particles  and  lumped  forces  acting 
on  the  particles,  the  formulation  is  written  in  a  vector  form.  One  of  the  distinct 


advantages  of  a  vector  formulation  is  that  the  code  development  is  considerably  simpler 
than  a  matrix  analysis. 

Note  that  the  equations  of  motion  for  each  mass  particle  are  essentially  solved 
independently.  Mass  particles  are  related  to  one  another  only  through  the  internal  forces. 
This  property  of  a  mass-force  system  offers  considerable  advantages: 

*  The  motion  of  each  particle  is  calculated  independently.  The  equations  of  motion  yield 
total  absolute  motion  of  the  particle,  including  both  the  rigid  body  motion  and  deformable 
motion. 

*  The  system  is  an  assemblage  of  independent  particles,  regardless  whether  these  particles 
form  a  unit  body  or  multiple  bodies.  The  only  difference  lies  in  the  calculation  of  internal 
forces.  Hence,  the  algorithm  automatically  permits  multiple  bodies  and  body  separations. 

*  Since  material  properties  are  included  in  the  computation  of  internal  forces,  inelastic  and 
discontinuous  material  properties  are  easy  to  handle. 

*  Adding  or  subtracting  a  set  of  nodes  do  not  affect  the  original  discretization.  Hence, 
creating  new  surfaces  or  eliminating  portions  of  the  medium  due  to  failure  presents  no 
numerical  problem.  Hence,  the  algorithm  can  be  modified  to  simulate  construction  process 
and  fracture. 

*  Numbering  of  the  nodes  has  no  effect  on  the  computation.  Again  this  is  convenient  for 
creating  or  eliminating  surfaces  and  bodies. 

*  Since  the  interactions  between  two  components  are  through  the  internal  forces,  the  co- 
existence of  soft  and  stiff  components  do  not  present  numerical  problems.  The  extreme 
case  is  that  some  components  have  zero  stiffness  as  a  result  of  failure. 

*  It  is  easy  to  handle  different  types  of  structural  component  without  the  need  for 
complicated  assemblage  process. 

*  Coding  for  the  algorithm  is  simple.  The  resulting  program  is  short  and  compact. 


Since  Newton's  law  of  motion  forms  the  basis,  the  algorithm  is  a  transient  procedure  by 
nature.  Static  solutions  can  be  obtained  by  attenuating  the  motion  through  the  use  of 
dynamic  relaxation  and  apply  external  loads  incrementally. 

b.  Explicit  Time  Integration: 

To  avoid  the  complexity  of  iterations,  a  simple  explicit  time  integration  formulation  is 
suggested  to  find  velocity  and  displacement  for  each  time  increment.  This  simplifies  the 
implementation  for  complicated  material  models,  changing  constraint  conditions  and 
loading  conditions.  Inelastic  and  failure  conditions  become  much  simpler  to  incorporate. 
However,  small  time  and  force  increments  are  required  for  numerical  stability.  This  leads 
to  longer  computational  times  in  general. 

c.  A  Co-rotational  for  Large  Rotations: 

To  develop  a  more  accurate  approach  to  handle  large  deformation  and  yet  simple  in 
computation,  a  co-rotational  formulation  of  the  kinematics  is  introduced.  For  each  time 
increment,  the  motion  is  separated  into  two  parts:  a  pure  rigid  body  rotation  of  the  entire 
finite  element  and  a  deformation  of  the  element  which  characterizes  the  shape  changes.  A 
convected  local  coordinate  system  is  introduced  which  is  updated  for  each  time  increment. 
The  strain  tensors  and  hence  the  corresponding  stress  tensors  are  formulated  in  the  local 
coordinates.  For  small  changes  in  the  element  geometry,  linear  stress  and  strain 
relationships  can  be  assumed  despite  the  large  overall  changes  in  geometry  due  to 
rotations. 


1.2  OBJECTIVES 

The  primary  objective  is  to  develop  a  three-dimensional  finite  element  program  for  the 
analysis  of  general  pavement  problems.  The  program  considers  conventional  static  and 
dynamic  loading  conditions  including  harmonic  excitations,  pulse  loadings,  ramp  loadings, 
and  multiple  step  loadings.  Provisions  are  made  for  the  convenience  of  handling  non- 
conventional  loadings  such  as  falling  weights  and  general  time-dependent  load  histories 
generated  from  non-destructive  testings  used  for  pavement  structural  evaluation.  One 
particular  important  aspect  of  the  program  is  to  develop  a  general  material  library. 
Common  soil  and  asphalt  material  models  in  the  form  of  linear  and  nonlinear  elastic 
materials,  elastic-plastic  materials  with  hardening,  and  viscoelastic  materials  are  included 
in  the  material  library. 

The  program  should  be  coded  such  that  further  expansions  of  the  load  and  material 
libraries  only  require  minimum  effort. 


1.3  LITERATURE  REVIEW 

a.  Algorithm  for  Pavement  analysis  and  design  using  3D  finite  elements: 

A  TRIS  data  base  search  was  performed.  Most  of  the  current  work  can  be  categorized  as 
the  use  of  finite  element  for  the  study  of  special  features  of  pavements,  rather  than  the 
development  of  a  suitable  finite  element  algorithm  for  the  pavement  analysis.  The  majority 
of  the  work  are  also  based  on  the  use  of  existing  commercial  codes. 


A  study  at  University  of  Minnesota  ( Koubaa  and  Krauthammer  1990)  was  reported.  The 
goal  is  to  develop  an  analysis  method  combined  with  a  non-destructive  testing  procedure 
for  the  evaluation  of  load  transfer  for  joints  in  concrete  pavements.  The  basic  analysis 
involved  a  frequency  response  analysis  by  dynamically  loading  the  joints.    A  three- 
dimensional  finite  element  method  was  used  to  analyze  various  joint  conditions  for  load 
transfer  ranging  from  full  to  partial  load  transfer. 

Stoner,  et  al.  1990  reported  on  research  plans  to  develop  a  3D  finite  element  program  to 
study  concrete  pavements  and  to  simulate  truck  actions.  The  objectives  were  to  develop  a 
truck  simulation  model,  to  develop  a  model  for  doweled  concrete  pavement,  and  to 
implement  these  two  models  in  an  interactive  fashion.  The  analysis  was  intended  to  obtain 
performance  relationships  based  on  damages  predicted  by  the  program.  Actual  damages 
recorded  on  Interstate  80  were  used  to  verify  the  predictions. 

Barksdale  1971  reported  a  study  of  compressive  stress  pulse  at  different  depths  in  a 
flexible  pavement.  Loadings  for  variable  vehicle  speed  were  considered.  A  series  of 
pseudo-dynamic  linear  and  nonlinear  elastic  analyses  were  conducted.  It  was  concluded 
that  linear  elastic  finite  element  was  adequate  for  the  analysis.  Dynamic  effects  including 
damping  and  inertial  forces  were  neglected  in  the  study  and  hence,  a  correction  factor  was 
introduced  in  order  to  match  the  results  of  AASHTO  road  tests. 

Paterson  1983  reported  on  the  finite  element  analysis  of  an  asphalt  overlay  of  a  cracked 
airport  concrete  pavement  in  combination  with  a  thin  interlayer  of  elastomeric  asphalt. 
Predictions  were  obtained  in  terms  of  an  equivalent  thickness  of  asphalt  overlay  which 
yields  the  same  performance  of  the  interlayer  system. 


Measured  pavement  deflections  in  combination  with  a  three-dimensional  finite  element 
analysis  were  used  to  evaluate  overlay  requirement  and  pavement  performance  in  a  study 
reported  by  Bala  and  Kennedy  1986.  The  reponse  predictions  of  the  finite  element 
analysis  were  calibrated  against  surface  deflections  measured  by  a  deflectograph.  The 
calibration  process  included  adjustments  determined  through  a  parametric  study.  The 
properties  assumed  for  the  materials  in  the  analysis  were  compared  with  laboratory 
determined  dynamic  test  results  for  in  situ  samples. 

Zaghloul  and  White  1993  a  reported  results  of  a  three-dimensional  dynamic  finite  element 
analysis  of  flexible  pavements.  The  analysis  simulated  actual  truck  loads  moving  at 
different  speeds.  Linear  and  non-linear  material  properties  were  used  to  model  different 
paving  materials  and  subgrades.  An  extended  Drucker-Prager  model  was  used  to  model 
granular  materials,  and  an  extended  Cam-Clay  model  was  used  for  the  clayey  soils. 
Asphalt  mixtures  were  modeled  as  viscoelastic  materials.  Including  these  material  models 
lead  to  the  capability  to  obtain  accurate  elastic  and  plastic  pavement  responses.  With  this 
capability,  they  are  able  to  predict  or  to  interpret  pavement  performances  under  a  variety 
of  loading  conditions  and  for  different  material  characteristics. 

The  3D  finite  element  analysis  was  verified  by  comparing  the  predictions  with  a  multi-layer 
elastic  system,  assuming  linear  elastic  properties  and  static  loads.  A  linear  correlation  was 
found  between  the  results  obtained  by  the  finite  element  predictions  and  the  multi-layer 
elastic  analysis.  To  verify  the  dynamic,  nonlinear  finite  element  analysis,  the  results  were 
compared  with  actual  measurements  of  pavement  deflection.  Agreement  at  95% 
confidence  level  was  obtained  between  the  deflection  predictions  and  the  measurements. 

A  sensitivity  study  was  performed  by  using  the  3D  finite  elements  for  the  effect  of  cross 
section  and  load  parameters  on  pavement  responses.  It  was  found  that  the  speed  of 


moving  vehicle  load  has  a  significant  effect  on  elastic  and  plastic  pavement  responses.  The 
confinement  of  shoulders  has  the  effect  of  reducing  pavement  deflections.  A  crack  along 
the  pavement/shoulder  joint  results  in  an  increase  in  pavement  deflection.  Temperature 
affects  the  asphalt  layer  and  hence  the  overall  pavement  responses.  The  loading  time  and 
the  rate  of  loading  were  found  to  have  significant  effect  also.  When  a  subgrade  is 
subjected  to  a  high  stress  level,  higher  than  its  yield  stress,  rutting  increases  significantly. 

The  effects  of  different  load  attributes,  axle  load  and  spacing,  number  of  axles,  and 
number  of  wheels,  as  well  as  cross  section  attributes,  subgrade  type,  material  properties 
and  the  type  of  deep  foundation  were  also  studied,  and  found  to  be  significant  to  pavement 
responses. 

In  a  separate  study  founded  by  INDOT,  Zaghloul  and  White  1993b  reported  the  study  of 
heavy  loads  and  their  effects  on  the  Indiana  highway  network.  It  was  decided  to  evaluate 
the  damages  induced  by  the  heavy  loads  by  their  LEF's.  Available  LEF's  were  found  not 
suitable  for  the  study.  New  LEF's  were  developed  based  on  the  total  pavement 
deformation  at  the  pavement  surface.  A  3D  dynamic  finite  element  model  was  used  to 
perform  the  analysis.  A  comparison  of  the  Purdue  LEF's  and  the  AASHTO  LEF's  was 
made  and  no  significant  difference  was  found.  Purdue  LEF's  consider  different  load  and 
cross  section  parameters,  while  the  AASHTO  LEF's  do  not.  In  addition,  Purdue  LEF's 
were  developed  based  on  analytical  models,  which  can  be  extended  to  consider  a  wide 
range  of  other  variables. 

The  study  reported  by  Zaghloul  and  White  have  utilized  a  commercial  general  purpose 
finite  element  code.  Due  the  general  nature  of  the  code,  the  run  time  required  to  use  a 
commercial  code  for  a  particular  problem  is  usually  much  longer  than  what  is  necessary. 


More  importantly,  the  flexibility  in  application  is  limited.  Modifications  to  tailor  the  code 
for  extended  pavement  applications  are  difficult. 

b.  Explicit  Finite  Element  Algorithms: 

Classical  static  solutions  of  finite  element  analysis  generally  adopt  an  implicit  approach. 
The  formulations  are  reduced  to  a  set  of  simultaneous  equations  (or  matrix  equations)  to 
find  the  displacement  and  stress  values  at  discrete  points.  The  details  of  an  implicit 
approach  are  well  documented  in  the  textbooks,  for  example,  (Bathe  1982).  To  extend 
the  formulations  to  handle  dynamic  analysis  due  to  impact  loadings,  or  to  perform 
nonlinear  analysis  due  to  inelastic  behaviors  of  material  and  large  changes  in  structural 
geometry,  one  encounters  considerable  difficulties  in  the  development  of  a  reliable  solution 
algorithm.  Existing  algorithms  generally  involve  complicated  integration  procedures  as 
well  as  repeated  iterations.  As  the  result,  the  application  of  an  implicit  approach  for 
general  transient  inelastic  analysis  is  somewhat  awkward.  Code  development  also  becomes 
difficult  and  inefficient  (Bathe,  et  al.  1975). 

To  circumvent  many  of  these  difficulties,  explicit  approaches  have  been  proposed  and 
implemented.  One  which  is  widely  used  considers  the  structure  as  a  set  of  discrete  mass 
particles.  By  using  the  particle  mass  equations  of  motion  including  external  and  internal 
forces,  the  formulation  reduces  the  finite  element  model  to  a  set  of  vector  equations. 
Thus,  the  approach  avoids  the  assemblage  and  the  storage  of  large  matrix  equations,  and 
thus  the  solution  algorithms  are  much  simplified.  The  vector  equations  generally  calculate 
particle  accelerations  due  to  external  loadings.  Displacements  and  corresponding  stresses 
can  be  found  by  a  time  integration  procedure. 
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The  explicit  approach  has  been  proved  to  be  particularly  suitable  for  impact  analysis  for 
which  a  time  series  of  the  loading  history  is  prescribed.  It  is  also  convenient  to  study 
highly  nonlinear  and  discontinuous  material  properties,  as  well  as  very  large  displacements 
often  associated  with  nonlinear  properties  (Key  1974,  Oden  and  Key  1973).  A  large 
number  of  general  purpose  computer  codes  have  been  developed  in  recent  years,  mostly 
tailored  for  applications  in  the  defense  industry  and  for  the  safety  analysis  of  nuclear 
power  plants.  For  example,  STRAW  was  developed  by  Argonne  National  Laboratory 
(Kennedy  et  al.  1985),  and  DYNA2D  and  DYNA3D  by  Lawrence  Livermore  Laboratory 
(Hallquist  1982).  Commercial  codes  including  ABAQUS  have  also  issued  explicit  versions 
recently. 

Fundamental  studies  of  explicit  approaches  have  been  carried  out  at  Purdue  Civil 
Engineering  during  the  past  decade.  Algorithms  for  large  deflection  of  space  frames  (Saha 
and  Ting,  1983),  for  concrete  fracture  analysis  (Saha  1983),  and  for  the  failure  analysis  of 
concrete  slabs,  folded  plates,  and  shells  (Labbane  1991)  have  been  reported.  They  include 
a  variety  of  structural  finite  element  formulations.  Special  algorithms  have  been  studied 
for  elastic- viscoplastic  materials  (Labbane  and  Ting  1991),  for  extremely  large  deflections 
(Rice  and  Ting  1991),  and  for  material  fragmentation  process  (Rice  and  Ting  1992). 
Computer  codes  using  the  explicit  approach  for  all  the  studies  have  shown  to  be  efficient 
and  short.  They  generally  require  a  small  memory  and  hence,  are  suitable  for  personal 
computers.  In  general  the  explicit  approaches  are  coded  for  transient  dynamic  analysis.  A 
dynamic  relaxation  procedure  is  used  to  find  static  solutions. 
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CHAPTER  2  CO-ROTATIONAL  APPROACH  AND  EQUATION  OF  MOTION 


In  this  chapter,  the  traditional  co-rotational  approach  for  a  plane  element 
subjected  to  displacements  is  formulated.  The  co-rotational  approach  is  used  to  treat 
large  rotations  of  a  deformable  solid  element. 

Let  the  solid  be  modeled  by  elements.  A  set  of  rectangular  coordinates  is 
attached  to  each  element,  which  translates  and  rotates  with  the  element  during  the 
deformation  process  but  does  not  deform.  Thus,  this  set  of  convected  or  co-rotational 
coordinates  is  always  an  orthogonal  system. 

If  the  solid  element  is  subjected  to  large  displacements,  we  assume  that  the 
translation  and  rotation  of  the  convected  coordinates  can  be  large.  However,  the 
deformations  described  in  the  convected  coordinate  system  are  small.  This  forms  a 
large  rotation-small  deformation  theory,  where  small  strain  theory  and  linear  stress- 
strain  relationships  can  be  used. 

Schematically,  we  may  illustrate  the  co-rotational  approach  by  considering  three 
c 
structural  configurations,  as  shown  in  Figure  2.1.     They  are  (a)  the  undeformed 

geometry  X  at  time  t  =  0,  (b)  the  deformed  geometry  x  at  time  r,  and  (c)  a  convected 

geometry  x  at  time  t .  Hence, 

dx  =  F  dX  ,      dx  =  Tdx 

T  and  F  are  transformation  matrices.  If  J  is  a  pure  rotation, 
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dx  =  TTdx 
where  TT  is  the  transpose  matrix  of  T  .  And 

dx  =  TF  dX  =  F  dX 

A 

The  strain  induced  by  the  transformation  F    can  be  treated  by  an  infinitesimal 
theory. 

Putting  in  the  context  of  finite  elements,  we  consider  a  nodal  displacement  vector 

a 

dt  in  global  coordinates  and  a  nodal  displacement  vector  d,.  in  convected  coordinates. 
They  are  related  by  a  rotation  matrix  T  as 

4  =  Td,  (2.1) 

We  further  assume  that  the  total  displacement  vector  can  be  written  as  the  sum  of  a 
deformation  vector  4   and  a  rigid  body  motion  vector  d/   , 

4  =  ded  +  dj  (2.2) 

The  corresponding  vector  in  convected  coordinates  are 

A  A(j  A  (. 

4  =  4+4  (2.3) 

and 

id  =  T4d  .  4r  =  T4 

For  small  deformation,  the  strain  tensor  depends  on  4  only. 
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2.1  Rotation  Matrix  for  Plane  Solid  Elements 


Consider  a  four-node  plane  element  shown  in  Figure  2.2(a).  (x,y)  is  a  set  of 
local,  or  convected,  coordinates  with  its  origin  at  the  center  of  the  element,  (x,  y)  is 
denoted  as  the  global  coordinates,  (u,  v)  is  the  displacement  vector  in  the  global 
coordinates.  From  time  t  to  time  t',  the  change  of  the  orientation  of  (x,y) 
coordinates  with  respect  to  the  rigid  body  rotation  angle  6  of  the  element  can  be 
described  as  shown  in  Figure  2.2(b)  such  as 


6  =— (a,  +a2  +a3  +a4) 


(2.4) 


or 


i    4 

e=-Ya(. 

4tr  ■ 


(2.5) 


where  a ,  is  the  rotational  angle  of  the  ith  local  nodal  position  vector  from  time  t  to 
time  t'.  By  using  cross  product,  a,  can  be  calculated  as 


Zx/'  = 


since,  k 


(2.6) 


or  a .  =  sin    < 


— [(*.-  -  xc  hi +vi-y'c)-  (y>  -  yc  X*.-  +  M<  -  *  e )] 


(2.7) 


where 


=  V(*.--*c)2+(y,--:>'c)2 


v 


=  ^(xi+ui-x'c)2+(yl+Vi-y'cf 


Xc  =-(^+X2+X3+Xi) 
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2.2  Kinematics  for  a  Plane  Solid 


The  displacement  vector  (  u,  v  )  can  be  decomposed  to  a  rigid  body  component 
and  a  deformation  component,  i.e., 

(2.13) 


u 

u' 

+ 

ud 

V 

v' 

V* 

Also,  the  rigid  body  vector  is  the  sum  of  a  translation  and  a  rotation  components, 


u' 

u' 

+ 

M° 

v' 

v' 

vy 

(2.14) 

as  shown  in  Figure  2.4.     Using  the  rotation  matrix,  R  ,  we  get  the  deformation 
component  in  convected  coordinates 


ud 

vd 


=    R 


u 


(2.15) 


assumed  an  infinitesimal  strain  theory. 

Note  that,  for  the  element  deformation,  the  rigid  body  motions  should  be 
separated  from  the  total  nodal  displacements.  Since  x  -axis  passes  through  the  element 

center,  the  nodal  displacement  components  da  and  dcy   of  the  central  point  of  the 

element  are  induced  by  rigid  body  motions.  This  deformation  displacement  vector 
(u  ,v  )  is  related  to  nodal  deformation  displacements  by  shape  functions.  In  finite 
element  technique,  the  isoparametric  formulation  then  is  introduced  for  this  term. 

Isoparametric  formulation  uses  the  same  shape  function  to  define  the  geometric 
shape  of  the  element  and  to  describe  the  displacements  within  the  element.  The 
function  is  formulated  using  a  natural  coordinate  system,  or  so-called  s-t  domain, 
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which  is  a  transformation  mapping  domain  of  the  convected  coordinates.  Since  there 
are  only  two  nodes  for  each  side  of  a  four-node  element,  a  linear  shape  function  for  the 
displacements  and  the  nodal  coordinates  along  the  s  or  t  element  boundary  is  assumed. 
Multiplication  of  the  linear  functions  in  s  and  t  yields  the  shape  functions  for  the  four- 
node  element  in  the  form 


u  =  a,  +  a2s  +  a3t  +  a4st 


v  =  a5  +  a6s  +  a7t  +  asst 


(2.16) 


x  =  a,  +  a2s  +  a3t  +  a4st 
y  =  a5+  a6s  +  ant  +  asst 


(2.17) 


By  solving  for  the  eight  a(  unknowns  in  terms  of  nodal  coordinates,  the  following 
equations  are  established 


*l- 

"*i 

0 

^2 

0 

^3 

0 

^4 

0  " 

y\~ 

0 

tfi 

0 

^2 

0 

^3 

0 

#4. 

X, 


y3 


>4j 


(2.18) 


l"l 
Ivl 


Nt      0      N2      0      Wj      0      N4      0 
0     N,      0     AT,      0     iV„      0      N* 


u. 


vu 


u. 


Lv4j 


(2.19) 
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A  schematic  diagram  for  the  transformation  mapping  is  shown  in  Figure  2.5,  and  the 
shape  function,  Nj ,  are  now 

Nx  =I(i-5)(i-r)       N2  =|(l+5)(l-0 

4  4 


N,  =±(l+s)(l+t)        N4  =  Ul-s)(l  +  t) 
4  4 


(2.20) 


which  satisfies  the  displacement  continuity  at  all  boundaries  and  Nj  +N2  +N3  +N4=l. 
The  element  strains  based  on  the  convected  geometry  are  given  by 


du 

£; 

dx 

X 

dv 

•  £- 

>  =  - 

■■■■■            S 

y 

^   . 

t*J 

3m    3v 

■■■       -j-  -  - 

dy     dx 

or  in  s-t  domain 


1 

3y  3m    3y  3m 
3r  35     35  3r 

0 
3Jc  3v     3r  3v 

0 

3ic  3v     dx  dv 

ds  3f     3f  35 
3y  3m     3y  3m 

_3j  3r     dt  ds 

3r  35     35  dt 

Substitute  Equation  (2.20)  into  Equation  (2.22b),  we  get 

e  =  Bd 

where  d  is  given  by  Equation  (2.10)  and 


(2.21) 


(2.22) 


(2.23) 


fi-r-r-jfij     B2     B3     BA 


(2.24) 


\J\  is  the  Jacobian  determinate  defined  as 
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Or  explicitly, 


u\= 


dx 

dy 

ds 

ds 

dx 

dx 

dt 

dt 

\A-U 


A|        X2        X^        X^ 


} 


0  l-t  t-s  s-l~ 

t-\  0  5+1  -S-t 

S-t  -5-1  0  t+l 

1-5  S  +  t  -t-\  0 


y~i 


(2.25) 


Also, 


5.  = 


dy  dNt 

dy  dNt 

o 

dt   ds 

ds   dt 
0 

d5t    dN; 

dx  dNt 

dxdNL 

dx  dNt 

ds   dt 
dy  dNt 

dt   ds 
dy  dNt 

ds   dt       dt   ds       dt   ds       ds   dt 


(2.26) 


and 


f^=  \[yl(s-l)  +  y2(.-l-s)  +  %(l  +  s)  +  y4a-s)] 

dt      4 


^  =  \[Mt-l)  +  y2(l-t)  +  y3(l  +  t)  +  y4(-l-t)] 

dx       1 

—  =  -[x1(t-l)  +  x2(l-t)  +  x3(l  +  t)  +  Xt(-l-t)] 
ds      4 

dx       1 

—  =  -[xl(s-l)  +  x2(-l-s)  +  x3(l  +  s)  +  x,a-s)] 
dt      4 


!H<-> 

3/       4 

lH<>-<> 

dr       4 

f=>> 

T3  4<1+*> 

3f       4 

^  =i(-'-o 

3"<4a-> 

ds       4 


8r       4 
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2.3  Principle  of  Virtual  Work 

Similar  to   the   traditional   finite  element   formulation,   the  equilibrium  of 
structural  system  is  defined  by  the  principle  of  virtual  work. 

l8Ue-X5W;  =  0  (2.27) 

where  the  summation  is  carried  out  for  all  the  elements  and  5C/eand  8W^are  the 
element  internal  and  external  virtual  work,  respectively. 

The  element  internal  virtual  work  in  terms  of  the  element  stress  and  strain  can 
be  written  in  the  form 

dUe=^SeT6  dV  (2.28) 

t  —      — 

For  infinitesimal  deformation, 

dV  =  dV0     and      Ve  =  V0 
where  V0  is  the  undeformed  volume.  Using 

S£T  =  (Sde)TBT 

Equation  (2.28)  becomes 

5Ue=(Sdef[   BT6dVQ  (2.29) 

0  ~ 

Since 

5Ue=$deffemt  (2.30) 

where  f™  is  the  internal  nodal  forces  associated  with  the  nodal  displacements. 
Hence, 
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femt=\BradVQ  (2.31) 

Equation  (2.31)  can  be  written  in  terms  of  the  material  stress-strain  relationships  which 
will  be  described  in  Chapter  four.  We  consider 

6=Ce  (2.32) 

where  C  is  the  material  property  tensor.  Then 

femt=(lBTCBdV0)de  (2.33) 

*»o    -    -  - 

For  a  two  dimensional  solid  with  thickness  t,  we  have 

dVQ  =  t  cIAq  =  t  dx  dy 
Equation  (2.33)  becomes 

feat=(t\JBjCBdxdy)de 

or  in  the  s- 1  domain, 

feM  =  (tflf_lBT  CB\j\ds  dt)de  (2.34) 

For  a  four-node  plane  element,  the  integration  of  Equadon  (2.34)  yields  eight  force 
components: 

f'r={L  K  L  fiy  L  hy  L  hy]       (2.35) 

Using  the  global  nodal,  displacements    de     through  the  transformation,  de  —Tde  , 

yields 

8  Ue  =  5  dj  TT  f*  =  8  dj  /,*"  (2.36) 
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where 


T    rint 


fr  =  t'  / 


(2-37) 


is  the  internal  nodal  force  vector  represented  in  the  global  coordinates. 

The  external  virtual  work  may  be  treated  as  the  sum  of  the  work  due  to  applied 
forces  and  the  work  due  to  inertia  forces.  For  the  inertia  forces,  if  the  element  mass  is 
lumped  as  concentrated  masses  located  at  the  nodes  of  the  element,  then  the 
corresponding  virtual  work  due  to  inertia  force  is 


WeA  =  -5  d[  Mede=S  d]  Me  de 


(2.38) 


where  Me  is  the  lump  mass  matrix  which  is  an  invariant  and  diagonal.  For  a  four-node 


plane  element, 
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0 

0 

0 

1 

(2.39) 


where  p  is  the  mass  density,  A  and  t  are  the  area  and  thickness,  respectively,  of  the 
element. 

The  virtual  work  due  to  applied  forces  can  be  formulated  following  the 
traditional  finite  element  approaches.  Briefly,  the  forces  given  in  global  coordinates 
should  be  transformed  to  components  in  convected  coordinates.    The  corresponding 
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nodal  forces  are  evaluated  according  to  the  shape  functions  along  the  element 
boundary.  After  sub-assembling  the  force  components  to  obtain  the  external  force 
vector  in  convected  coordinates,  it  is  transformed  back  to  global  coordinates  and  we 
get 

W.'-*dIfr  (2.40) 

where  //"   is  the  external  nodal  force  vector  in  global  coordinates  for  the  element. 

To  consider  the  body  force,  since  the  element  mass  is  lumped  in  the  nodes,  it  is 
convenient  to  treat  the  body  force  as  a  concentrated  external  nodal  force  with  respect 
to  the  gravity  and  to  be  included  in  Equation  (2.40)  as  part  of  the  //*  . 

The  total  external  virtual  work  for  the  element  is  then 

SWe=$WeA+5WeB  =8dJ(f<«-Medt)  (2.41) 

Using  the  principle  of  virtual  work,  we  get 

Z5Ue-X5\^  =  0 


or 


2ZdIfr-l5dJ(fr-Mede)  =  0  (2.42) 

e         -       -         e  -      - 

Introducing  a  global  (or  assembled)  nodal  displacement  matrix  d,  Equation  (2.42) 
becomes 

§dT{F'M-Faa+Md)  =  0  (2.43) 
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where  FmI  ,  Fext  and  M  are  assembled  internal  nodal  force  matrix,  external  force 
matrix,  and  mass  matrix,  respectively.  Or 

e 

F^^/r  (2.44) 

e 
e 

The  square  matrix  M  remains  to  be  diagonal  after  assemblage.  As  each  of  5  d  is 

arbitrary,  we  get 

FBj"-FB"  +  M«da  =  0  ,    a=l,2,3, n 

where  n  is  the  total  number  of  unknown  nodal  displacements.  Or,  the  equadons  of 
modon  are 

■  da  =  j^  (  FaCT  -  Faint )  (2.45) 

Note  that  due  to  the  use  of  lumped  masses,  the  element  formulation  only  requires  the 
assemblage  of  force  vectors.  Furthermore,  the  equation  of  motion  for  each  unknown 
can  be  calculated  individually  without  the  need  for  solving  simultaneous  equation. 
Thus,  the  current  formulation  combined  with  an  explicit  time  integration  to  solve  the 
equations  of  motion  presents  significant  advantages  in  code  development,  data  storage, 
and  overall  efficiency. 
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-*X 


v 

(a) 

Figure  2.1  Three  Structural  Configurations  in  Co^rotational 
Approach  (a)  the  Undeformed  Geometry  at  time  t  =  0,  (b)  the 
Deformed  Geometry  at  time  t,  (c)  a  Convected  Geometry  at 
time  t. 
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-»»  X 

(a) 

Figure  2.2  (a)  Displacement  and  Rotation  of  Lines  of  a  Plane 
Element  in  the  X-Y  Plane,  (b)  Pure  Rotation. 
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Figure  2.3  Degrees  of  Freedom  of  Nodes  of  a  Four-Node 
Plane  Element  in  X-Y  Plane  :  (a)  Displacement  Part,  (b)  Force 
Part. 
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Figure  2.4  Decomposition  of  Displacements  (a)  Translation  and 
Rotation  without  Strain,  (b)  Pure  Strain  without  Translation  and 
Rotation. 


27 


Figure  2.5  (a)  Square  Element  in  s-t  Coordinates,  (b)  Quadrilateral 
Element  in  X-Y  Coordinates  are  Mapped  into  a  Square  Element  in  s-t 
Coordinates. 
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CHAPTER  3  EXPLICIT  TIME  INTEGRATION 


Direct  integration  procedures  commonly  used  in  solving  structural  dynamic  problems 
may  be  categorized  as  explicit  or  implicit  methods.  The  main  difference  is  that  explicit 
methods,  such  as  central  difference  method,  calculate  the  displacements  at  time  t  +  At 
based  on  the  equation  of  motion  at  time  t ;  while  the  implicit  methods,  such  as  Houbolt, 
Wilson-9  and  Newmark-{3  method  [Bathe,  1982],  use  the  equation  of  motion  at  time  t  + 
At  .  Hence,  iterations  are  generally  required  for  implicit  procedures. 

To  choose  a  suitable  method,  one  should  consider  the  combined  effect  of 
discretizations  on  time  and  space,  for  example,  calculating  the  natural  frequencies  of  a 
structure.  In  the  spatial  discretization,  if  the  mass  matrix  is  obtained  by  a  lumped 
approach,  the  frequencies  are  often  underestimated.  If  however,  the  mass  matrix  is  found 
by  assuming  a  consistent  approach,  the  frequencies  are  often  overestimated.  Similarly,  in 
the  time  discretization,  an  explicit  integration  causes  the  frequencies  to  be  overestimated; 
while  the  implicit  integration  underestimates  the  frequencies  [Key,  1978].  Therefore,  it 
seems  that  a  favorable  combination  is  to  choose  an  explicit  time  integration  method  with 
lumped  masses,  or  an  implicit  method  with  consistent  masses. 

An  advantage  of  using  an  explicit  method  with  lumped  matrices  is  that  there  is  no 
need  to  assemble  stiffness  matrices  and  does  not  need  to  solve  simultaneous  equations. 
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The  solutions  can  be  carried  out  in  vector  form  which  requires  very  small  storage  space. 
The  drawback  is  that  the  procedure  is  conditionally  stable.  Hence,  a  small  time  increment 
is  required  for  the  calculation. 

For  a  multi-degree-freedom  system,  the  time  integration  technique  used  in  this  work 
is  a  second  order  central  difference  formulation  which  has  the  following  relationships  for 
accelerations  and  velocities  : 

I  =-iI(da+1-2dn  +  dn_1)  (3.1) 

At 
1 

4   =   2^(4n-<k-l)  (3-2) 

where  At  is  the  time  increment,  and  the  time  span  t  =  n  At  . 
The  equation  of  motion  at  the  n  -  th  time  step  has  the  form 

d,  =  M-'CFr-F,"")  (3.3a) 

If  we  wish  to  find  a  quasi-static  solution  through  a  dynamic  relaxation  procedure,  a 
damping  force  may  be  added, 

d„  =  M-'CF^-Fr-F^)  (3.3b) 

The  damping  force  may  be  written  by  assuming  a  standard  Rayleigh  damping  [James, 
1994] . 

Fndmp=  Cd,  =  (<xM  +  PK)4 

in  which  a  and  (3  are  constants,  K  is  the  stiffness  matrix.  In  our  calculation,  since  the 
global  stiffness  matrix  is  not  available  in  the  formulation,  (3  is  assumed  to  be  zero.  Or, 
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Fadmp  =  Ci  =  ccMi  (3.4) 

Substituting  Equation  (3.1),  (3.2)  and  (3.4)  into  (3.3b)  yields 

^  =(_2_)^l(Fr  _ir-)  +  (_i_yn  -(IZ^^  ,  (3.5) 

al       2  +  aAt    M      "  2  +  ocA/     "      2  +  aA/    "  ' 

Note  that  using  Equation  (3.5),  the  displacements  at  t  +  At  are  calculated  by  using  the 
mass  values  and  the  external  and  internal  forces  of  the  previous  time  step.  An  important 
simplification  can  be  introduced  by  assuming  a  diagonal  mass  matrix,  since  the  inverse 
matrix  can  be  obtained  by  the  reciprocals  of  the  diagonal  mass  values.  All  the  calculations 
in  Equation  (3.5)  involve  vector  operations  only. 

To  start  the  solution  process  to  calculate  d,  ,  we  need  d_,  .  Using  the  initial 
condition  dg  ,  we  may  write 

d0  =  ^(d1-d.I) 
Or, 

d_,  =  d,-2Atdo  (3.6) 

Substituting  Equation  (3.6)  into  Equation  (3.5),  d,  can  be  solved  for. 
Alternatively,  we  may  use  Equation  (3.2)  which  gives 

d,  =  2At  do  +  d_,  (3.7) 

Substituting  Equation  (3.7)  into  Equation  (3.1)  for  n  =  0  gives 


^  =  -L(dI-2d0  +  d.I) 

At 


Or, 
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22 
d_,  =  do  -At  do  +  ^At24  (3.8) 

where  d^  can  be  found  from  Equation  (3.3b),  such  as 

4  =  M-,(F0--F,h-aM4)  (3.9) 

Since  the  time  increment  is  small  relative  to  the  time  span  of  interest,  there  is  no  significant 
difference  in  these  two  starting  procedures. 

In  general,  the  explicit  time  integradon  is  conditionally  stable.  That  is,  the  time 
increment  has  to  be  smaller  than  a  limit  to  avoid  calculations  becoming  divergent.  For  a 
multi-degrees  of  freedom  system,  such  a  limit  is  difficult  to  obtain.  Hence  ,  an 
approximate  limit  on  At  is  suggested  [Hughes,  1979] 

^  <  -^sr  (3.10) 

with  CO  "■*  being  the  highest  frequency  value  for  the  elements  in  a  structure. 
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CHAPTER  4     MATERIAL  MODELS  AND  STRESS-STRAIN  RELATIONSHIPS 

In  this  Chapter,  the  basic  material  models  which  have  been  implemented  and  verified  in  the 
current  computer  codes  are  described.  The  stress-strain  relationships  are  briefly  reviewed. 
They  include: 

a.  Linear  Elastic  Models: 

Stress-strain  relationships  are  given  in  the  form  of  Hooke's  Law.  The  elastic  constants  are 
Young's  modulus  and  Poisson's  ratio. 

A  generalization  is  to  input  tangent  modulus  to  replace  the  Young's  modulus.  Effectively, 
this  is  a  bi-linear  elastic  model  or  a  nonlinear  elastic  model. 

b.  Plastic  Models: 

Stress-strain  relationships  are  given  for  the  Drucker-Prager  yield  criterion  with  associated 
flow  rules.  Three  hardening  rules  are  implemented:  isotropic,  kinematical,  and  mixed 
hardening  rules.  The  special  case  of  selecting  a  Mises  yield  criterion  for  the  yield 
criterion  and  assuming  an  associated  flow  rule  is  included.  Another  alternative  which  is 
often  used  to  model  soil  behaviors  is  to  select  a  Drucker-Prager  yielding  function  and 
assuming  a  non-associated  flow  rule  with  Mises  function  as  the  plastic  potential.  This 
combination  has  also  been  implemented. 

Material  constants  inputs  for  soils  are  specified.  They  are  the  cohesion  and  the  internal 
frictional  angle. 

c.  Viscoelastic  Models: 

Stress-strain  relationships  are  given  for  the  viscoelastic  material  of  Maxwell  type. 
Implementation  for  other  types  of  viscoelastic  model  such  as  the  standard  linear  solid  and 
Burger  type  are  trivial.  A  general  creep  model  based  on  a  curve-fitting  of  the  creep  data 
can  be  included  by  a  simple  modification  of  the  viscoelastic  material  subroutine. 
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4.1  Linear  Elastic  Stress-Strain  Relations 

A  material  possesses  a  linearly  elastic  property  if  it  is  elastic  and  the  one-to-one 
relationship  between  the  stresses  and  strains  is  linear.  The  generalized  Hooke's  law 
expressing  the  stress  components  as  linear  homogeneous  functions  of  the  strain 
components  is  given  as 

*q  =  CijldZu  (4.1) 

where  CijU  is  the  fourth  order  material  property  tensor  and  can  be  expressed  in  terms 

of  two  independent  material  properties  if  material  is  isotropic,  which  means  properties 
are  not  changed  by  orthogonal  transformations  of  coordinates. 

For  isotropic  case,  the  material  property  tensor  can  be  expressed  as 

Cm  =  *8  #8  u  +  p.  (8  tt8  n  +  5  „5  jk )  (4.2) 

where  8  ..  is  the  Kronecker  delta,  X  and  ll  are  Lame  constants  in  terms  of  Young's 
modulus  E  and  Poisson's  ratiov  or  shear  modulus  G  as 

x  _  Ev  E 

(l+v)(l-2v)         ^  2(l+v)  (4-3) 

Plane  stress  case  assumes  C  3J  =  0,  whereas,  plane  strain  case  assumes  £  3  =  0. 

In  geotechnical  field,  plane  strain  condition  is  more  common. 

4.2  Plastic  Stress-Strain  Relations 

Plasticity  theory  extends  elasticity  theory  when  the  state  of  stress  satisfies  the 
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yield  criterion.  It  offers  a  mathematical  description  of  the  mechanical  behavior  of 
material  in  the  plastic  range.  Since  soils  behave  dependent  on  volumetric  stress,  the 
pressure-dependent  plastic  model,  including  Mohr-Coulomb  and  Drucker-Pracker  yield 
criteria,  are  more  popular  for  granular  soils.  In  this  research,  the  Drucker-Prager 
criterion  is  chosen.  It  is  mentioned  that  in  this  section  the  sign  convention  for  stresses 
is  defined  as  "positive  means  compression"  following  the  soil  mechanics  convention. 


4.2.1  Drucker-Prager  Yield  Criterion 

The  Drucker-Prager  yield  function  is  given  as  the  follows 

/(/[  72,a,K)  =  J2  -(a/j  +K)2  =0  (4.4) 

where  I{  =  the  first  invariant  of  stress  tensor  =  G  .. , 

J 2  =  the  second  invariant  of  stress  de viator  tensor  =— G  °  G  f  ,  in  which 

°ij  =«J//-j/i6(/,and 

OC  ,K  =  material  constants. 

In  terms  of  the  more  familiar  soil  parameters,  cohesion  c   and  friction  angle  <{>  ,  the 

material  constants  a  and  K  can  be  taken  either  one  set  of  the  follows 

g  _         2sin(j)  6c-cos<j) 

V3  (3  +  sin  <|> )  ~  VJ (3  +  sin  <j) )  (4'5) 

n             2sin(l>              „               6c-cos6 
or,  a  -  —7=— and  K  = " 


V3(3-sin<|))  V3(3-sin<|>)  (4"6) 
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Figure  4. 1  shows  the  yield  surface  in  the  Haigh-Westergaard  stress  space  and  on  the  Pi- 
plane,  respectively. 


4.2.2  Incremental  Stress-Strain  Relationships 

In  plastic  region,  Equation  (4.1)  is  no  longer  valid  and  is  revised  as 

c  a  =  c'/ki  £  u  (4.7) 

where  C'fkl  is  the  elastic-plastic  material  property  tensor  and  can  be  calculated  from 
the  elastic  material  property  tensor,  Cljkl ,  and  the  plastic  material  property  tensor, 
cm  .  by 

CW  ~  Cw  ~  cuki  (4.8) 

Note  that  the  plasdc  stress  tensor  is  function  of  the  current  stresses,  a  .. ,  independent 

of  the  increment  of  stresses,  da  i}  .  To  find  C£, ,  the  basic  ingredients  of  the  theory  of 
plasticity  must  be  carried  out  based  on  the  given  a  tj  ,  Cijkl ,  dzkl,  and/. 

5/, 


/da 

(1)  Unit  normal:  n-  = 


^     PA 


'da  , 


(4.9) 


i.         \df  /     II      r  df      df      V  a. 

Wh£re    I  A  it  (da-dt)     =  the  ma^de  of  df/l 


.  /  dGa 


If  stresses 


undergo  plastic,  two  conditions  must  be  satisfied:  (a)  /  =  0  and  (b)da?  nr>0t 


36 


where  dc'^  =  Cijkld£ld 


(2)  Strain  decomposition: 


de  ■■  =  de '-  +  d e  I 
"*w  i]  ij  tj 


(4.10) 


in  which  the  incremental  small  strain,  dz  .■ ,  is  decomposed  into  an  elastic  part,  dz  ?• , 


and  a  plastic  part,  de  » 


(3)  Constitutive  equation: 


(4)  Flow  rule: 


da  ij  =  Cijki  dfL  li 


da/ 

OCT   ;; 


/OCT.,. 


(4.11) 


(4.12) 


where  a    is  the  consistency  parameter  which  is  a  scalar,  and  q  is  the  plastic  potential 
function.    If  q  =  f ,  this  is  so-called  associated  flow  rule  condition;  if  q  &  f  ,  non- 


associated  flow  rule. 


(4.13) 


(5)  Consistency  condition:  /=0  and  df=0. 

A  schematic  consistency  condition  is  shown  in  Figure  4.2. 

Note  that  the  use  of  the  associated  flow  rule  with  Drucker-Prager  model  often 
overestimates  the  plastic  dilation  of  the  soil  (shown  in  Figure  4.3);  therefore,  it  is  thus 
common  to  use  a  non-associated  flow  rule  with  Drucker-Prager  to  correct  the 
problem.  The  von  Mises  yield  function  is  then  given  as  the  plastic  potential  while 
Drucker-Prager  yield  function  is  for  the  yield  surface.  The  von  Mises  yield  function 
can  be  generated  from  the  Drucker-Prager  yield  function  by  setting  c  £  0  and  fy  =  0 . 
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4.2.2.1  Linear  Elastic-Perfectly  Plastic  Model 

If  the  stress  strain  relationship  of  the  material  is  modeled  as  an  elastic-perfectly 
plastic  behavior,  shown  in  Figure  4.4,  the  elastic -plastic  material  property  tensor  for 
Drucker-Prager  yield  function  with  non-associated  flow  rule,  as  described  previously, 
can  be  derived  based  on  Equation  (4.9)  to  (4.13)  and  concludes  the  forms 


Drucker-Prager  function:  /  (Ix  72,CC,k)  =  J2  -  (cc/j  +k)    =0 


(4.14) 


von  Mises  function:      q(J2  ,K)  =  J2  —  K     =0 


A  = 


n  :  C  :  dz 

n  :  C  :  m 


(4.15) 
(4.16) 


yda       C 


m  = 


h. 


dc 


•V2K 


5/ 


da       aD-2a(cc/1+K)l 


n  = 


5/ 


do 


aD-2a(a/[+K)l 


(4.17) 


(4.18) 


da  =  Cep:dz  =  (C-  Cp):dz 


(4.19) 


C:m,®  ^:  C 

c  =  ^-= ^^ 

-  n  :  C  :  m 


(4.20) 


where  1  is  the  second  order  unit  tensor  with  components  L  =  5 , 


!/  >J 
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4.3  Viscoelastic  Stress-Strain  Relations 

Viscoelastic  models  characterize  material  behaviors  which  are  time-dependent  and 
temperature-dependent  Uniaxial  tests  or  simple  shear  tests  are  commonly  used  to  rind  the 
material  constants.  Hence,  viscoelastic  models  are  generally  determined  for  uniaxial  (Young's 
modulus)  or  shear  (shear  modulus)  stress-strain  relationships. 

Extending  the  models  for  bi-axial  or  tri-axial  stress  conditions,  assumptions  have  to  be 
made.  Two  cases  which  are  commonly  assumed  for  application  are  taking  the  material  to  be 
incompressible  and  assuming  the  Poisson's  ratio  to  be  a  constant 

In  the  program,  we  have  taken  Poisson's  ratio  to  remain  constant  and  the  Young's 
modulus  is  replaced  by  a  viscoelastic  model. 

4.3.1   Viscoelastic  Models 

Material  characterized  by  two  constants  are  known  to  be  a  Maxwell  type,  the  relaxation 
function  (Young's  modulus)  has  the  form 


E-Ee-'  <4'21> 

o 


where  x  is  a  relaxation  time. 

If  four  material  constants  are  used,  it  is  known  as  a  Burger's  model,  where 


E  =  E,  e  "*  ♦  E,  e  *  (4'22) 


4.3.2  Creep  Models 

Using  the  creep  test  data,  time  dependent  behavior  can  be  characterized  by  creep  models 


formulated  in  incremental  form.   Let  incremental  creep  strain  Aec  be 
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A£_  = 


'<JV 


=  fo4n 

where  n,  m,  k  are  constants  obtained  from  curve-fitting. 
Then,  the  incrementation  stress  Aa  is 

Aa  =  E(Ae  -  Aec) 

Ae 
=  E(l-—L)Ae 
Ae 


(4.23) 


(4.24) 


Implementation  of  the  creep  models  is  trivial,  merely  modify  the  tangent  modulus  E  for  each 
time  or  load  increment. 
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Figure    4.1     Drucker-Prager    Yield    Surface 
Westergaaxd  Stress  Space,  (b)  on  71  Plane. 


(a) 


in    Haigh- 
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Figure  4.2  Schematic  Consistency  Condition 


Figure  4.3  The  Dilation  of  Drucker-Prager  Yield  Function  with 
Associated  Flow  Rule 
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Figure  4.4  Idealized  Stress-Strain  Curves:         Elastic-Perfectly 
Plastic    Model 
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CHAPTER  5.   THREE  DIMENSIONAL  SOLID  ELEMENTS 

In  chapter  2,  general  formulations  of  the  co-rotational  approach,  the  principle  of  virtual 
work  and  the  general  equations  of  motion  are  described.  Specific  formulations  for  a  plane  solid 
element  are  given. 

For  a  three-dimensional  solid  element,  the  formulations  for  the  rotation,  virtual  work,  and 
the  equations  are  the  same.  The  only  differences  are: 

(a)  there  are  three  rotational  matrices  for  the  transformation  between  the  convected 
coordinates  and  the  global  coordinates, 

(b)  there  are  three  displacement  components  and  six  independent  stress  and  strain 
components. 

5.1  Rotational  Matrix  for  3D  Solid  Elements 

The  transformation  between  a  set  of  convected  coordinates     (x,yji)     and  the  global 

coordinates  (x,y,z)  requires  three  angles,  known  as  the  Euler  angles  in  calculus.  Computation 
of  the  angles  is  the  same  as  given  in  Eq.  (2.5).  For  each  plane,  (xy),  (y,z)  and  (zx),  a  rotational 
angle  can  be  calculated  which  yields  Euler  angles  (9Z,  8X,  6y).  Then,  the  rotation  matrix  for  the 
transformation  is 


*=*x*y*z  (5-1) 


where    R      R      and  R        are  rotation  matrix  for  a  plane  solid. 
-x  ,~y  ,        ~z 
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For  example, 


R 


c    -s 

0 

s    c 

0 

0    0 

1 

(5.2) 


where  c=  cos  0X  and  s=  sin  0X 

5.2  Isoparametric  Elements  for  3D  Solid 

An  eight-node  3D  isoparametric  element  is  considered.  There  are  three  nodal 
displacements  (u^Vj.Wj)  at  each  node.  Hence,  the  element  has  24  degrees  of  freedom.  The  shape 
functions  are 


x  =  ^Nfatrfx 


i=l 


~  I 


(5.3) 


«  =  5>.(s,r,r)« 


;=i 


(5.4) 


where  x  = 


and  u  = 


(s,  t,  r)  are  three  natural  coordinates.  The  eight  shape  functions  have  the  form 


N.  =  i.  (1  +  ss)  (1   +  tt)  (1  +  rr)        i  =  1,2,3,...,8 

o 


(5.5) 


where  (%,  tj,  r^  are  position  values  of  the  node  in  the  natural  coordinates.    The  six  strain 


components  can  be  calculated  by 
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a/* 

written  in  the  natural  coordinates,  the  derivate    _    for  example,  becomes 

dx 


dx         \J\ 


fs  y, 

*■> 

ft  y, 

z« 

fr  y, 

2, 

(5.7) 


where  fs  =    _  ,  etc. 
ds 


The  Jacobian  IJI  is 
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171  = 


x,  y,  *t 

xr  yr  z, 


(5.8) 


Fig.  5.1  shows  the  sketch  of  an  eight-node  isoparametric  element. 
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Figure  5.1     An  Eight-Node  Three-Dimensional  Isoparametric  Solid  Element 
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Conclusions 

1.  The  algorithm  which  is  based  on  a  vector  formulation  of  finite  element  for  the 
equations  of  motion,  an  explicit  time  integration  technique  for  the  solutions,  and  a  co- 
rotational  formulation  for  the  large  displacement  seem  to  work  well.  The  programs  are 
efficient  for  transient  dynamic  problems.  However,  static  solutions  are  also  shown  to  be 
reliable,  and  accurate. 

2.  The  programs  based  on  this  algorithm  are  compact,  efficient,  and  flexible  for 
expansions. 

3.  A  well-known  general  purpose  commercial  code  ANSYS  is  used  to  verify  the  codes. 

4.  Six  verification  problems  are  considered  for  the  three  dimensional  code  and  five  for  the 
two  dimensional  code. 

5.  Material  models  included  in  the  codes  are  elastic,  elastic-plastic  (three  versions),  and 
viscoelastic.  Modifications  of  the  material  library  are  convenient. 

6.  The  load  library  has  included  general  time  history,  ground  acceleration,  and  harmonic 
inputs.  By  specifying  the  time  histories,  pulse  loadings  and  ramp  loadings  for  static 
solutions  can  be  considered. 
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Appendix  1 

•  A  summary  of  the  programs 

•  Using  Microsoft  Excel  to  prepare  input  data  and  to  plot  output  data 

•  User's  Manual  for  a  two  dimensional  finite  element  program  Solid2D  (S2DP) 

•  User's  Manual  for  a  three  dimensional  finite  element  program  Soild3D  (S3DP) 

•  Fortran  Source  Code  for  Solid2D 

•  Fortran  Source  Code  for  Solid3D 
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A  SUMMARY  OF  THE  PROGRAMS:  Solid2D  and  Solid3D 

1.  Load  library 

The  current  versions  have  essentially  programmed  for  three  types  of  forces  applied  at  the 
nodes. 

a.  An  arbitrary  time  history  of  the  force 

b.  A  ground  acceleration  input 

c.  A  sinusoidal  force  with  fixed  frequency  and  amplitude. 

Currently,  it  is  set  up  for  6  different  types  of  force  input  or  6  types  of  acceleration  input. 
They  can  be  applied  at  30  different  nodes.    For  the  sinusoidal  force,  it  can  be  applied  at  10 
nodes.  The  user  may  input  500  time  increments  for  the  time  histories  of  acceleration  and 
forces.  These  numbers  can  easily  be  increased  by  changing  the  dimension  statements  in 
one  subroutine  "readata". 

Pulse  impact  forces  and  variable  step  loadings  are  prescribed  by  the  force  values  at 
different  time  increments.  To  obtain  static  solutions,  force  input  is  carried  out  by 
prescribing  load  history  as  an  one-step  time  history,  a  multi-step  load  history  or  a  ramp 
load. 

2.  Boundary  and  Initial  conditions 

Displacement  constraints  are  prescribed  at  the  nodes. 

By  default,  the  initial  velocities  the  initial  displacements  are  zero.  That  is,  the  structure  is 
at  rest  initially.  The  programs  allow  the  initial  displacements  and  initial  velocities  to  be 
prescribed.  Currently,  they  may  be  prescribed  at  10  different  nodes.  To  increase  the 
number,  merely  change  the  dimension  sizes  in  the  subroutine  "readata". 
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3.  Body  Forces 

Gravity  acceleration  can  be  input. 

4.  Material  Library 

To  retain  maximum  flexibility  for  the  future  development,  the  material  data  inputs  are  not 
classified  according  to  each  individual  material  model.  Thus,  material  models  can  easily  be 
combined  to  represent  emerging  developments  in  material  modeling.  For  each  material 
type,  all  the  material  data  can  be  made  available. 

Current  material  inputs  are:  Mass  density,  Young's  modulus,  Poisson's  ratio, 

Tensile  strength  (uniaxial  tensile  yield  stress), 

Cohesion  and  Internal  angle  of  friction  (for  soil  only), 

Tangent  modulus  (slope  of  the  stress  and  strain  curve  after  yielding), 

Hardening  rule  (assuming  linear  work-hardening),  and 

Relaxation  time  (assuming  Maxwell  model). 

The  material  models  available  are: 

1 .  linear  elastic  material 

2.  linear  viscoelastic  material  of  Maxwell  type 

3 .  elastic-plastic  material  with  associated  flow  rule  assuming  Mises  Criterion 

4.  elastic-plastic  material  with  associated  flow  rule  assuming  Drucker-Prager  Criterion 

5.  elastic-plastic  material  with  non-associated  flow  rule 
(Drucker-Prager  for  yielding  and  Mises  for  hardening) 

For  all  the  elastic-plastic  materials,  there  is  a  choice  of  kinematical  hardening  rule, 
isotropic  hardening  rule,  and  mixed  hardening  rule  (characterized  by  a  factor  beta 
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as  the  percentage  of  kinematical  rule). 

With  minor  modifications,  bi-linear  elastic  material,  nonlinear  elastic  material,  complex 
viscoelastic  models,  nonlinear  creep  models,  and  viscoelastic-plastic  material  are  readily  to 
be  implemented. 
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Using  Microsoft  Excel  to  prepare  input  data  and  to  plot  output  data 


How  to  create  the  Input  Data  File 

In  case  of  small  and  simple  finite  element  mesh,  DOS  editor  is  preferred  to  use  to  create  input 
data  file.  However,  if  a  larger  or  complicate  mesh,  the  MS  Excel  is  more  preferable  that  the  previous  one. 
In  the  next  following  section,  the  procedure  of  creating  input  data  file  by  MS  Excel  is  illustrated. 

1.  Follow  the  Input  Guide,  and  start  from  card  1  (each  data  number  occupying  one  cell  in  MS  Excel) 

2.  When  you  get  to  card  4(coordinate  of  each  node),  you  can  just  enter  information  for  the  first  node  and 
use  copy  command  to  do  the  other  nodes 

3.  Following  the  same  procedure  for  generating  the  element  mesh(card  5). 

4.  For  card  6  to  card  21,  enter  the  value  of  each  data. 

5.  Save  the  input  data  file  in  the  DOS  format  under  's2dp.dat  —  for  solid2d  and  s3dp.dat  —  for  solid3d'. 

How  to  Plot  the  Results  in  the  MS  Excel 

1.  Open  the  output  file  or  's2dp.out  or  s3dp.out'  in  the  MS  Excel. 

2.  By  selecting  file  type  as  'fixed  width'  and  click  'next' 

3.  Select  'column  break  line'  for  each  particular  data  set  and  click  'next'. 

4.  Select  'data  column  format'  as  'general'  and  then  click  'finish' 

5.  Select  the  sets  of  data  to  be  plotted  by  shading  the  interested  area. 

6.  Click  Chart  Wizard  and  locate  your  plotted  area. 

7.  Select  the  chart  type(recommend  the  XY  scatter  or  Line). 

8.  Select  the  format  for  the  chart 

9.  Following  the  command  on  the  screen  to  select  the  chart  set  up. 

10.  Click  Finish. 
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Solid2D  User  Manual 

and 

Input  Data  Guide 


************************************************************************ 


INPUT  GUIDE (for  Solid2D  program) 
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*  Sub.  dynamic  : 

*  1.  head  (20a4) 


head    :  problem  description  (  80  characters  ) 
2 .  iprob, nnd, nel , nummat , numout , ndof , maxstp, delta , alpha 


iprob 

nnd 

nel 

nummat 

numout 

ndof 

maxstp 

delta 

alpha 

toler 

gravity: 


no.  of  time-steps  skipped  between  outputs 

total  no.  of  nodal  points 

total  no.  of  elements 

total  no.  of  different  materials  or  sections  used 

total  no.  of  output  records  requested  each  d.o.f. 

of  each  node's  d,v,a,sigma  forms  each  output  record. 

degree  of  freedom  per  node  (=2  ) 

max.  time  steps  or  cycles  of  calculation 

time  increment  (sec) ,  must  be  less  than  critical  value 

coeff.  of  mass  damping 

tolerance  limit  as  a  switch  to  stop  running  program 

when  the  increment  of  displacement  is  smaller  than  it 

acceleration  of  gravity 

(1)  no  gravity  load  imposed  if  gravity=0 . 

(2)  Gravity  load  is  imposed  in  the  initial  condition. 

(3)  Displacement  due  to  the  gravity  load  imposed  is 
NOT  initialized. 


3 .  iacc, iforce, inital, imesh, iplane 


* 
* 

iacc 

* 

* 

iforce 

* 

* 

* 

* 

* 

inital 

* 

* 

* 

imesh 

* 

* 

* 

iplane 

* 

index  for  ground  acceleration  time-history  input 

0  =  no  ace.      1  =  x-dir      2  =  y-dir 
index  for  force  function 

0  =  no  external  force  function 

1  =  arbitrary  shape  force  function 

2  =  type  of  F=  f*sin(wt) 

3  =  type  of  F=  f*cos(wt) 
index  for  initial  condition 

0  :  d=0  v=0        2  :  d=0  v 

1  :  d    v=0        3  :  d    v 

index  for  printing  out  nodal  coor.  &  element  data 

0  =  not  print  out 

1  =      print  out 
index  for  plane  stress  =  1 

or  plane  strain  =  2 


*  Sub, 

*  4. 
* 

* 

* 

* 
* 

* 


readata  : 
n,xc (n) ,yc (n) , (ifixx(i) , i=l, 2) 


n 


node  no. 


xc(n)  :  x 
yc(n)  :  y 
ifixx(l) :  x 
if ixx(2) :  y 


coord,  of  node  n 
coord,  of  node  n 
translation  B.C. 
translation  B.C. 


(  0  =  free 


fixed  ) 


5 .  n , node ( 1 ) , node ( 2 ) ,  , node ( 7 ) , node ( 8 ) , knt 
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n 
node (1) 
node ( 5 ) 
node ( 6 ) 
node (7) 
node (8) 

knt 


node ( 4 )  :  node  4 


element  no. 

node  1  for  element  n    ==: 

material  property  no. 

row  no.     for  element  n 

column  no.  for  element  n 

element  condition  (=1)  (for  checking  Jacobian  use  ) 

no.  of  node  number  increment  used  for  automatic 

element  generation  (for  horizontal  &  vertical  dir. ) 


6.  k,e(l,k),e(2,k) e(ll,k) 


k 
e(l,k) 
e(2,k) 
e(3,k) 
e(4,k) 
e(5,k) 
e(6,k) 


material  group  no. 

material  type  no . 

mass  density 

young ' s  modulus 

poisson  ratio 

tensile  strength  of  metal 

cohesion  for  soil 


(  = 


(  = 


1 — metal  ;=  . 
rho  ) 

eyng  =  E  ) 
poisson   ) 
ft  ) 
cohesion  ) 


soil) 


7.  e(7,k)  ,e(8,k)  , ,e(12,k) 


e(7,k)  :  friction  angle  for  soil 
e(8,k)  :  material  model 


e(9,k)   :  tangent  modulus 
e(10,k)  :  hardening  rules 

(  =  0 
e(ll,k)  :  plate  thickness 
e(12,k)  :  tau  for  Viscoelastic 


(  skip  8,9,10  ,if  iacc  =  0  ) 
nacc,npnts 


(  =  phiangle  ) 
(  =  0  for  Linear  Elastic 
=  1  for  Viscoelastic 
=  2  for  von  Mises, 
=  3  for  Drucker-Prager) 
(  =  eynt  =  Et  ) 
(  =  beta   )  (0  =<  beta  =<  1) 
kinematic  ;  =  1  :  isotropic  ) 


9. 

10. 


nacc   :   total  no.  of  different  ground  acceleration  history 
npnts  :   no.  of  time-acc.  pairs  in  each  ace.  history 

g     :   gravity  acceleration  (unit  must  be  consistent) 

ta ( j , i) ,aa( j , i)  j=l, npnts   i=l,nacc 


ta(j,i)  :  time   for  ace. 
aa(j,i)  :  value  for  ace. 


11. 


(  skip  11  —  16  ,if  iforce  =  0  )  -- 
(  skip  11  --  14  ,if  iforce  .ne.  1  ) 
numif ,nnaf 


12. 


13. 


numif  :   total  no.  of  impact  force  history 

nnaf   :   total  no.  of  nodes  applied  by  arbitrary  shape  impact 
force  function 

kfpnts(i)  :  total  no.  of  time-force  pairs  in  a  force  function 
history   [  1<  kfpnts(i)  =<  500  ] 

tf ( j , i) , ff ( j , i)  j=l,jf   i=l, numif   jf=kfpnts(i) 

tf(j,i)  :  time   for  force 
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*  ff(j,i)  :  value  for  force  * 

*  * 

*  14.  jnode(i) , jdir (i) , jaf (i)  i=l,nnaf  * 

*     -  * 

*  jnode(i)  :  node  no.  where  impact  force  is  applied  * 

*  jdir(i)  :  d.o.f.  corresponding  to  which  this  force  is  applied  * 

*  1  =  x-dir  2  =  y-dir  * 

*  jaf (i)  :  time-force  history  no.  * 

*  * 

*  15.  nnaf   :  total  no.  of  nodes  applied  by  sinusoidal  impact  force  * 

*     function  * 

*  * 

*  16.  jnode (i) , jdir (i) , f (i) , omega (i)  i=i,nnaf  * 

*  * 

*  jnode (i) :  node  no.  where  impact  force  is  applied  * 

*  jdir(i)  :  d.o.f.  corresponding  to  which  this  force  is  applied  * 

*  1  =  x-dir  2  =  y-dir  * 

*  f(i)  :  amplitude  of  this  sinusoidal  force  F=  f*sin(wt)  * 

*  omega (i):  frequency  "       "              F=  f*cos(wt)  * 

*  * 

*  (  skip  17  --  20  ,if  inital  =  0  )  * 

*  17.  ndisi  :  total  no.  of  displacement  type  I.e.  * 

*     * 

*  18.  ndnod(i) , kdis (i) ,disi (i)  i=l, ndisi  * 

*     * 

*  ndnod(i) :  node  no.  where  displ  type  I.e.  is  applied  * 

*  kdis(i):  d.o.f.  corresponding  to  which  this  I.e.  is  applied  * 

*  1  =  x-dir  2  =  y-dir  * 

*  disi(i):  value  of  initial  displacement  * 

*  * 

*  19.  nveli  :  total  no.  of  velocity  type  I.e.  * 

*     * 

*  20.  nvnod(i) , kvel (i) ,veli (i)  i=l, nveli  * 

*     * 

*  nvnod(i) :  node  no.  where  velicity  type  I.e.  is  applied 

*  kvel(i):  d.o.f.  corresponding  to  which  this  I.e.  is  applied  * 

*  1  =  x-dir  2  =  y-dir  * 

*  veli(i):  value  of  initial  velocity  * 

*  * 

*  21.  kout (1) ,kout (2) ,kout (3)  do  loop  i=l,numout  * 

*      * 

*  kout(l)  :  node  no.  for  which  response  output  is  requested  * 

*  kout (2)  :   =  0    for  displacement  request  * 

*  =1    for  velocity  request  * 

*  =2    for  acceleration  request  * 

*  =3    for  stresses  request  * 

*  kout (3)  :  global  dof  corresponding  to  which  output  is  requested  * 

*  1  =  x-dir  2  =  y-dir     3  =  xy-plane  * 

*  (for  shear  stress  only)  * 
************************************************************************ 

*  Limitation  * 

*  * 

*  1.  ta(500, 6) ,aa (500, 6)  only  for  6-type  ground  accelerations  * 

*  &  500  time-acc.  pairs  for  each  type  * 

*  2.  tf (500,6) ,ff  (500,6) only  for  6-type  force  functions  * 

*  &  500  time-force  pairs  for  each  type  * 

*  3.  jnode (30) , jdir (30) , jaf (30)  only  for  30  positions  applied  by  * 

*  each  type  force  function  * 

*  4.  jnode (30) , jdir (30) , f (10) , omega (10)  only  for  10  positions  * 

*  applied  by  sinusoidal  force  function  * 
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*  5.  ndnod(lO) ,kdis (10) ,disi (10)  only  for  10  positions  applied    * 

*  by  initial  displacement  * 

*  6.  nvnod(lO) , kvel (10) , veli (10)  only  for  10  positions  applied     * 

*  by  initial  velocity       -  * 

*  * 

*  If  the  limitations  were  violated  ,  please  change  the  dimension  of   * 

*  above  variables  in  sub."  readata  ".  * 

*  * 
************************************************************************ 
************************************************************************ 

*  Note  * 

*  * 

*  There  are  three  files  in  this  program:  * 

*  1.  Input  data  file  ====>     's2dp.dat'  * 

*  2.  Print  out  of  input  file     ====>     's2dp.in'  * 

*  3.  Output (result)  file  ====>  's2dp.out'  * 
************************************************************************ 
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SolicBDUser  Manual 

and 

Input  Data  Guide 
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*  *  *  *  * 

* 
* 


******************************************************************* 


INPUT  GUIDE (for  SOlid3D) 


*  Sub 

*  1. 
* 

* 

* 

*  2. 

* 

* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 

*  3. 


.  dynamic  : 
head  (20a4) 


head 


problem  description  (  80  characters  ) 


iprob , nnd, nel , nummat , numout , ndof , maxstp, delta , alpha 


total  no. 
total  no. 


no.  of  time-steps  skipped  between  outputs 
total  no.  of  nodal  points 

of  elements 

of  different  materials  or  sections  used 

total  no.  of  output  records  requested  each  d.o.f. 

of  each  node's  d,v,a,sigma  forms  each  output  record. 

degree  of  freedom  per  node  (=3  ) 

max.  time  steps  or  cycles  of  calculation 

time  increment  (sec) ,  must  be  less  than  critical  value 

coeff.  of  mass  damping 

tolerance  limit  as  a  switch  to  stop  running  program 

when  the  increment  of  displacement  is  smaller  than  it 

acceleration  of  gravity 

(1)  no  gravity  load  imposed  if  gravity=0. 

(2)  Gravity  load  is  imposed  in  the  initial  condition. 

(3)  Displacement  due  to  the  gravity  load  imposed  is 
NOT  initialized. 


iacc, iforce, inital, sideL, imesh, iplane 


iprob 

nnd 
nel 

nummat 
numout 

ndof 

maxstp 

delta 

alpha 

toler 

gravity: 


iacc    :  index  for  ground  acceleration  time-history  input 

0  =  no  ace.      1  =  x-dir   2  =  y-dir   3  =  z-dir 
iforce  :  index  for  force  function 

0  =  no  external  force  function 

1  =  arbitrary  shape  force  function 

2  =  type  of  F=  f*sin(wt) 

3  =  type  of  F=  f*cos(wt) 
inital  :  index  for  initial  condition 

0  :  d=0  v=0       2  :  d=0  v 

1  :  d   v=0       3  :  d   v 

imesh   :  index  for  printing  out  nodal  coor.  &  element  data 

0  =   not  print  out 

1  =      print  out 


*  Sub. 

*  4. 

* 

* 

* 
* 
* 
* 

* 
* 
* 


readata  : 
n,xc (n) ,yc (n) , zc (n) , (ifixx(i) , i=l, 3) 


n 
xc  (n) 
yc(n) 
zc  (n) 

ifixx(l) 
if ixx(2 ) 
if ixx (2 ) 


node  no. 

x  -  coord,  of  node  n 
y  -  coord,  of  node  n 
z  -  coord,  of  node  n 
x  -  translation  B.C. 
y  -  translation  B.C. 
z  -  translation  B.C. 


(  0  =  free 


1  =  fixed  ) 


5.  n,node (1) ,node (2) ,  , node (11) ,node (12) , knt 
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n 
node ( 1 ) 
node ( 9 ) 
node (10) 
node (11) 
node (12) 

lent 


node ( 8 )  :  node  8 


element  no. 

node  1  for  element  n   ==: 

material  property  no . 

row  no.     for  element  n 

column  no.  for  element  n 

element  condition  (=1)  (for  checking  Jacobian  use  ) 
no.  of  node  number  increment  used  for  automatic 
element  generation  (for  horizontal  &  vertical  dir.) 


6.  k,e(l,k)  ,e(2,k)  , ,e(6,k) 


k 
e(l,k) 
e(2,k) 
e(3,k) 
e(4,k) 
e(5,k) 
e(6,k) 


material  group  no . 

material  type  no. 

mass  density 

young ' s  modulus 

poisson  ratio 

tensile  strength  of  metal 

cohesion  for  soil 


7.  e(7,k),e(8,k) e(10,k) 


* 

* 
* 

*  if 

*  if 
* 


e(7,k)  :  friction  angle  for  soil 
e(8,k)  :  Material  Model 


(e(8,k)  =  0)  >  e(9..11,k)  =  0 

(e(8,k)  =  1)  >  e(9..10,k)  =  0 

e(9,k)  :  tangent  modulus 
e(10,k)  :  hardening  rules 

(  =  0 
e(ll,k)  :  Tau  for  Viscoelastic 


(=  l--metal;=  2 

(  =  rho  ) 

(  =  eyng  =  E  ) 

(  =  poisson   ) 

(  =  ft  ) 

(  =  cohesion  ) 


soil) 


(  =  phiangle  ) 

(  =  0  for  Linear  Elastic 
=  1  for  Viscolastic 
=  2  for  von  Mises, 
=  3  for  Drucker-Prager) 


(  =  eynt  =  Et  ) 
(  =  beta   )  (0  =<  beta  =<  1 
kinematic  ;  =  1  :  isotropic 


(  skip  8,9,10  ,if  iacc  =  0  ) 
nacc,npnts 


nacc   :   total  no.  of  different  ground  acceleration  history 
npnts  :   no.  of  time-acc .  pairs  in  each  ace.  history 

9.    g     :   gravity  acceleration  (unit  must  be  consistent) 

10.  ta ( j , i) ,aa ( j , i)  j=l, npnts   i=l,nacc 


ta(j,i)  :  time   for  ace. 
aa(j,i)  :  value  for  ace. 


11. 


(  skip  11  - 
numif ,nnaf 


16  ,if  iforce  =  0  ) 


total  no.  of  impact  force  history- 
total  no.  of  nodes  applied  by  arbitrary  shape  impact 
force  function 

:  total  no.  of  time-force  pairs  in  a  force  function 
history   [  1<  kfpnts(i)  =<  500  ] 


13.  tf ( j , i) , f f  ( j , i) j=l,jf   i=l, numif   jf=kfpnts(i) 


numif 
nnaf 


12.  kfpnts(i) 


tf(j,i)  :  time   for  force 
ff(j,i)  :  value  for  force 
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*  14.  jnode(i) , jdir(i) , jaf (i)  i=l,nnaf 


*  jnode(i)  :  node  no.  where  impact  force  is  applied  * 

*  jdir(i)  :  d.o.f.  corresponding  to  which  this  force  is  applied  * 

*  1  =  x-dir     2  =  y-dir    3  =  z-dir  * 

*  jaf (i)  :  time-force  history  no.  * 

*  * 

*  15.  nnaf   :  total  no.  of  nodes  applied  by  sinusoidal  impact  force  * 

*     function  * 

*  * 

*  16.  jnode (i) , jdir (i) , f (i) , omega (i)  i=i,nnaf  * 

*     * 

*  jnode (i):  node  no.  where  impact  force  is  applied  * 

*  jdir(i)  :  d.o.f.  corresponding  to  which  this  force  is  applied  * 

*  1  =  x-dir     2  =  y-dir    3  =  z-dir  * 

*  f(i)  :  amplitude  of  this  sinusoidal  force  F=  f*sin(wt)  * 

*  omega (i) :  frequency     "      "              F=  f*cos(wt)  * 

*  * 

*  (  skip  17  --  20  ,if  inital  =  0  )  * 

*  17.  ndisi  :  total  no.  of  displacement  type  I.e.  * 


*     * 

*  18.  ndnod(i) ,kdis (i) , disi (i)  i=l, ndisi  * 

*     * 

*  ndnod(i) :  node  no.  where  displ  type  I.e.  is  applied  * 

*  kdis(i):  d.o.f.  corresponding  to  which  this  I.C.  is  applied     * 

*  1    =   x-dir             2    =  y-dir           3    =    z-dir  * 

*  disi(i):  value  of  initial  displacement  * 

*  * 

*  19.  nveli  :  total  no.  of  velocity  type  I.C.  * 

*     * 

*  20.  nvnod(i) , kvel (i) , veli (i)  i=l;nveli  * 

*     * 


*  nvnod(i):  node  no.  where  velicity  type  I.C.  is  applied  * 

*  kvel(i):  d.o.f.  corresponding  to  which  this  I.C.  is  applied 

*  1  =  x-dir     2  =  y-dir    3  =  z-dir 

*  veli(i):  value  of  initial  velocity 


* 


* 


* 


* 
* 

*  * 

*  21.  kout (1) ,kout (2) ,kout (3)  do  loop  i=l,numout 

*     * 

*  kout(l)  :  node  no.  for  which  response  output  is  requested  * 

*  kout (2)  :   =  0    for  displacement  request  * 

*  =1    for  velocity  request 

*  =2    for  acceleration  request  * 

*  =3    for  stresses  request 

*  kout (3)  :  global  dof  corresponding  to  which  output  is  requested  * 

*  1   =  x-dir  2    =  y-dir            3    =   z-dir  * 

*  4  =  xy-plane   5  =  yz-plane   6  =  zx-plane  * 
•A********************************************************************** 

*  Limitation  * 

*  * 

*  1.  ta(500, 6) ,aa(500, 6)  only  for  6-type  ground  accelerations  * 

*  &  500  time-acc.  pairs  for  each  type  * 

*  2.  tf (500, 6) , ff (500, 6)  only  for  6-type  force  functions  * 

*  &  500  time- force  pairs  for  each  type  * 

*  3.  jnode (30) , jdir (30) , jaf (30)  only  for  30  positions  applied  by  * 

*  each  type  force  function  * 

*  4.  jnode(30) , jdir (30) , f (10) ,omega(10)  only  for  10  positions  * 

*  applied  by  sinusoidal  force  function  * 

*  5.  ndnod(10) ,kdis (10) ,disi (10)  only  for  10  positions  applied  * 

*  by  initial  displacement  * 
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*  6.  nvnod(lO) ,kvel (10) , veli (10)  only  for  10  positions  applied    * 

*  by  initial  velocity  * 

*  * 

*  If  the  limitations  were  violated  ,  please  change  the  dimension  of    * 

*  above  variables  in  sub."  readata  ".  * 

*  * 
************************************************************************ 
************************************************************************ 

*  Note  * 

*  * 

*  There  are  three  files  in  this  program:  * 

*  1.  Input  data  file  ====>     's3dp.dat'  * 

*  2.  Print  out  of  input  file      ====>      's3dp.in'  * 

*  3.  Output (result)  file         ====>     's3dp.out'  * 

************************************************************************ 
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Solid2D  Source  Code 
(S2DP.FOR) 
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************************************************************************ 

*  * 

*  Program  for  Dynamic  Plane  Solid  Problem  (  solid2D  )  * 

*  * 
************************************************************************ 

*  * 

*  This  program  is  for  analyzing  a  2-D  dynamic  problem  and  developed   * 

*  on  the  basis  of  :  * 

*  (1)  Traditionally  Co-Rotational  Approach  * 

*  (2)  Explicit  Time  Integration  Method  (central  difference)  * 

*  (3)  Lumped  Mass  Modeling  * 

*  (4)  4-node  Solid  Isoparametric  Element  * 

*  (5)  4 -point  Gaussian  Integration  * 

*  (6)  Elastic-lenear  work-hardening  Model  * 

*  (7)  Viscoelastic  Model  * 

*  (8)  von  Mises  yield  criterion  * 

*  (9)  Drucker-Prager  yield  criterion  * 

*  * 

*  Tat sana  Nil award  * 

*  Purdue  University  * 

*  03/07/96  * 

*  * 

*  * 
************************************************************** 

*  * 

*  Program  for  Dynamic  Plane  Solid  Problem  * 

*  (Solid2d  =  S2DP  )  * 

************************************************************************ 

c 

program  solid2D 

implicit  real*8  (a-h,o-z) 

dimension  ar (30000) 

maxq  =  30000 
c   create  the  input  and  output  filenames  (  extension  part  of  filename  -- 
c   — .dat  &  .out  ,will  be  created  automatically  ) 
c 

open  (5,  file= ' s2dp.dat ' ) 

open  (6,  f ile= ' s2dp. in' ) 

open  (7,  file= ' s2dp.out ' ) 

call  dynamic  (ar,maxq) 
c 

print  *,   'TOTALLY  COMPLETE  !  ' 

close  (unit=5) 

close  (unit=6) 

close  (unit=7) 
c 

stop 

end 
************************************************************************ 

*  * 

*  set  index  no.  of  each  variable  by  dynamic  allocation  method 

*  read  in  control  parameters   &   call  main  subroutines  * 

*  * 

************************************************************************ 

c 

subroutine  dynamic  (ar,maxq) 
c 

implicit  real*8  (a-h,o-z) 
dimension  ar (maxq) 
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dimens  i  on  head (20) 
c 

common  /box  1/  iprob, delta, alpha, toler, gravity 

common  /box  5/  maxcyc,maxout,maxmat 
c 

do  100  i=l,maxq 
100  ar(i)  =  0.0 
c 

read (5, 130)   head 

write (6, 140)  head 

write(6,150) 
c 

c   iprob  is  number  of  time-steps  skipped  between  output 

c 

read ( 5 , * )     iprob, nnd, nel , nummat , numout , ndof , maxstp , delta , alpha 

*  , toler, gravity 
c 

write (6,160)  iprob, nnd, nel, nummat, numout, ndof , maxstp, delta, alpha 

*  , toler, gravity 
c 

read (5, *)     iacc, i force, inita, imesh, iplane 

write (6, 170)  iacc, i force, inital, imesh, iplane 
c 

c   dynamic  allocation  

c 

meq    =  nnd*ndof 

maxdof  =  meq 

maxel   =  nel 

maxmat  =  nummat 

maxnod  =  nnd 

npint   =  1 

nxmass  =  npint+meq 

na      =  nxmass +meq 

nv     =  na+meq 

nd     =  nv+meq 

nxc    =  nd+meq 

nyc    =  nxc+maxnod 

nforce  =  nyc+maxnod 

nifix   =  nforce+meq 
c 

c   element  node  connectivity 
c 

nnode   =  nifix+meq 
c 

c  material  properties 
c 

nem    =  nnode+8*maxel 
c 

c   output  record 
c 

maxcyc  =  maxstp 

maxout  =  numout 
c 

if  (maxstp  .eq.  0)  maxcyc  =  1 

if  (numout  .eq.  0)  maxout  =  1 
c 

nkout   =  nem  +12*maxmat 

ndp     =  nkout  +  3*maxout 

ndn    =  ndp  +  meq 
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c 
c 


c 
c 


c 

c 
c 


nsigmaP  =  ndn  +  meq 
nsigmaN  =  nsigmaP  +  4*6*nel 
nepslonP  =  nsigmaN  +  4*6*nel 
nepslonN  =  nepslonP  +  4*6*nel 
nPLalphaP  =  nepslonN  +  4*6*nel 
nPLalphaN  =  nPLalphaP  +  4*6*nel 
nPLrP  =  nPLalphaN  +  4*6*nel 
nPLrN  =  nPLrP  +  4*nel 
nelplas  =  nPLrN  +  4*nel 
nelpout  =  nelplas  +  4*nel 
nsl  =  nelpout  +  4*nel 
ns2  =  nsl  +  nnd 
ns3  =  ns2  +  nnd 
ns4  =  ns3  +  nnd 
nproutv  =  ns4  +  nnd 
maxindex  =  nproutv  +  maxout 

i  f  (  maxindex  . gt .  maxq  )  then 
print  *,  '  There  is  not  enough  dimension  available.  ' 
print  *,  '   ==>  Please  increase  no.  in  ar ( . . . )  &  maxq= 

stop 

endif 


call  readata  (nnd, nel, nummat, numout, iacc,ndof, ar (nnode) , iforce, 
+  imesh, ar (nxc) , ar (nyc) , ar (nif ix) , ar (nem) , ar (nkout) , 

+  inital) 


print  *,  'call  readata  --  complete!' 


call  bmass  (nel,nnd,ndof ,ar (nem) ,ar (nnode) ,ar (nxc) ,ar (nyc) 
y  ,ar (nxmass) ) 


print  *,  'call  bmass 
time=0 . 0 


complete! ' 


call  esolv  (time, nel, nnd, ndof ,  iaccnumout,  iforce,  ar  (nxmass)  , 

+  ar(nforce) ,ar(npint) ,ar(nifix) ,ar(nd) ,ar(nv) ,ar(na) , 

+  ar(nxc) ,ar(nyc) , ar (nnode) ,ar(nem) , ar (nkout) ,maxstp, 

+  ar (ndn) , ar (ndp) , ini tal, meq, ar (nsigmaP) , ar (nsigmaN) , 

+  ar (nepslonP) , ar (nepslonN) , ar (nPLalphaP) , ar (nPLalphaN) , 

+  ar(nPLrP) ,ar(nPLrN) , ar (nelplas) , ar (nelpout) , ar (nproutv) 

+  ar (nsl) , ar (ns2) , ar (ns3) ,ar (ns4) , iplane) 


c 
c 


print 


130  format 
140  format 
150  format 
160  format 

+ 
+ 
+ 
+ 


'call  esolv 


complete ! 


(20a4) 

(/2x, 'card  l',5x,20a4) 

(lx,80( '-' ) ) 

(2x, 'card  2 ', 5x, 'parameter  card',/, 
15x, 'no  of  time-steps  skipped  between  outputs 
15x, 'number  of  nodes         =',il0,/, 
15x, 'number  of  elements      =',il0,/, 
15x, 'number  of  materials     =',il0,/, 


,i6, /, 
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,il0,/, 
,110,/, 
,il0,/f 

,el0.3,/, 
,el0.3,/, 
,el0.3,/, 
,fl0.5,/) 


+  15x, 'number  of  output  req 

+  15x, 'no.  of  d.o.f/node 

+  15x, 'no.  of  time  steps 

+  15x, 'time  increment 

+  15x, 'coeff  of  mass  damping 

+  15x, ' tolerance  limit 

+  15x, 'acceleration  of  gravity 
170  format  (2x, 'card  3',5x, 'index  card',/, 

+  15x, ' index  for  accel .       =',il0,/, 

+  15x, 'index  for  force       =',il0,/, 

+  15x, 'index  for  I.  C.       =',il0,/, 

+  15x, ' index  for  mesh  output (1)  or  not(0)       =',14, /, 

+  15x, ' index  for  plane  stress(l)  or  strain(2)    =',i4) 
c 

return 

end 

c 
a*********************************************************************** 

*  * 

*  calaulate  translational  mass  of  all  nodes  in  x-dir  &  y-dir  * 

*  and  form  Mass-matrix  all  at  one  time  * 

*  * 
a*********************************************************************** 

c 

subroutine  bmass  (nel,nnd,ndof , e, rnode,xc,yc,xmass) 
c 

implicit  real*8  (a-h,o-z) 

dimension  mode  (8,1)  ,  xc  ( 1 )  ,  yc  (1)  ,  e  ( 12  , 1 )  ,  xmass  ( 1 ) 

dimension  xmasc(8) 

common  /Gauss/  s (4) , t (4) ,w(4) 
c 

c   calculate  mass  for  each  node 
c 

meq  =  nnd*ndof 
c 

do  100  i=l,meq 
100  xmass(i)  =  0.0 
c 

c   set  Gauss  points  and  weights  (  w=l  for  4  point  Gauss  Quadrature  ) 
c 

const  =  dsqrt(3.0d0) 

s(l)  =  -1.0 /const 

t(l)  =  -1.0 /const 

w(l)  =   1.0 

s(2)  =   1.0 /const 

t(2)  =  -1.0/const 

w(2)  =   1.0 

s(3)  =   1.0/const 

t(3)  =   1.0/const 

w(3)  =   1.0 

s(4)  =  -1.0/const 

t(4)  =   1.0/const 

w(4)  =   1.0 
c 

c   loop  through  all  element 
c 

do  200  i=l,nel 
c 

c   set  node  no.  and  node  coord,  for  element 
c 
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c 
c 
c 


nl  =  int  (mode  (1,  i)  ) 

n2  =  int  (mode  (2,  i)  ) 

n3  =  int  (mode  (3,  i)  ) 

n4  =  int  (mode  (4,  i)  ) 

xl  =  xc(nl) 

x2  =  xc(n2) 

x3  =  xc(n3) 

x4  =  xc(n4) 

yl  =  yc(nl) 

y2  =  yc(n2) 

y3  =  yc(n3) 

y4  =  yc(n4) 

material  properties  for  element  (  b=thickness,  rho=mass  density  ) 

mtyp  =  int  (mode  (5,  i)  ) 
rho  =  e  ( 2  ,  mtyp ) 
b  =  e  ( 1 1 , mtyp ) 


c 
c 
c 


c 
c 
c 


do  300  j=l,8 
300  xmasc(j)  =0.0 

set  shape  function  Ni 

do  400  k=l,4 
shpfl  =  (1.0-s(k) ) 
shpf2  =  (1.0+s(k)) 
shpf3  =  (1.0+s(k) ) 
shpf4  =  (1.0-s(k) ) 


derivatives  of  shape  function  w.r.t.  s,t 


(1 

0-t(k)) 

/ 

4 

0 

(1 

0-t(k)) 

/ 

4 

0 

(1 

0+t(k)) 

/ 

4 

0 

(1 

0+t(k)) 

/ 

4 

0 

slds  = 

-d 

0-t(k)) 

/ 

4 

0 

s2ds  = 

(1 

0-t(k)) 

/ 

4 

0 

s3ds  = 

(1 

0+t(k)) 

/ 

4 

0 

s4ds  = 

-d 

0+t(k)) 

/ 

4 

0 

sldt  = 

-d 

0-s(k)) 

/ 

4 

0 

s2dt  = 

-d 

0+s(k) ) 

/ 

4 

0 

s3dt  = 

(1 

0+s(k) ) 

/ 

4 

0 

s4dt  = 

(1 

0-s(k)) 

/ 

4 

0 

c 
c 

derivatives 

of  global 

cood.  x,y  w.r.t.  s,t 

c 

xs  =  slds 

*xl  +  s2ds*x2  +  s3ds*x3  +  s4ds*x4 

ys  =  slds 

*yl  +  s2ds*y2  +  s3ds*y3  +  s4ds*y4 

xt  =  sldt 

*xl  +  s2dt*x2  +  s3dt*x3  +  s4dt*x4 

yt  =  sldt 

*yl  +  s2dt*y2  +  s3dt*y3  +  s4dt*y4 

c 

detJ  =  xs 

*yt  -  ys*xt 

c 
c 

compute  translational 

masses  xmasc (8) 

c 

xmasc (1) 

=  xmasc (1 

+  b*rho*w(k) *shpfl*detJ 

xmasc (2) 

=  xmasc (2 , 

+  b*rho*w(k) *shpfl*detJ 

xmasc (3 ) 

=  xmasc (3 , 

+  b*rho*w(k) *shpf2*detJ 

xmasc (4) 

=  xmasc (4) 

+  b*rho*w(k) *shpf2*detJ 

xmasc (5) 

=  xmasc (5 

+  b*rho*w(k) *shpf3*detJ 

xmasc (6) 

=  xmasc (6( 

+  b*rho*w(k) *shpf3*detJ 

xmasc (7) 

=  xmasc (71 

+  b*rho*w(k) *shpf4*detJ 

xmasc (8) 

=  xmasc (8' 

+  b*rho*w(k)*shpf4*detJ 
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c 
c 
c 


c 
c 
c 


400  continue 


calculate  index  of  mass  storage 


ml 
m2 
m3 
m4 


(nl-l)*ndof 
(n2-l)*ndof 
(n3-l)*ndof 
(n4-l)*ndof 
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xmass (ml+1) 
xmass (ml+2 ) 
xmass (m2+l) 
xmass (m2+2) 
xmass (m3+l) 
xmass (m3+2) 
xmass (m4+l) 
xmass (m4+2) 
continue 


1  mass  matrix 

xmass (8) 

xmass (ml+1 

+ 

xmasc ( 1 ) 

xmass (ml+2 

+ 

xmasc (2) 

xmass (m2+l 

+ 

xmasc (3) 

xmass (m2+2 

+ 

xmasc (4) 

xmass (m3+l 

+ 

xmasc (5) 

xmass (m3+2 

+ 

xmasc ( 6 ) 

xmass (m4+l 

+ 

xmasc (7) 

xmass (m4+2 

+ 

xmasc ( 8 ) 

** 

* 

* 

* 

* 

** 

c 


return 

end 
********************************************************************** 

* 
compute  elastic  constitutive  coefficients  under   3D  condition       * 

* 
********************************************************************** 

subroutine  elastd3d  (eyng,poisson,ed) 

implicit  real*8  (a-h,o-z) 
dimension  ed(6,6) 

do  50  i=l,6 

do  60  j=l,6 
ed(i,j)  =  O.dO 
60  continue 
50  continue 

G   =  eyng/ (2 .d0* (l.dO+poisson) ) 

dla   =  poisson*eyng/ ( (1 .d0+poisson) * (l.d0-2 .d0*poisson) ) 

+  dla 


ed(l 

1) 

=  2.d0- 

"G 

ed(2 

2) 

=  ed(l 

1) 

ed(3 

3) 

=  ed(l 

1) 

ed(4 

4) 

=  G 

ed(5 

5) 

=  ed(4 

4) 

ed(6 

6) 

=  ed(4 

4) 

ed(l 

2) 

=  dla 

ed(2 

1) 

=  ed(l 

2) 

ed(3 

1) 

=  ed(l 

2) 

ed(3 

2) 

=  ed(l 

2) 

ed(l 

3) 

=  ed(l 

2) 

ed(2 

3) 

=  ed(l 

2) 

return 

end 
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c 

************************************************************************ 

*  * 

*  Solve  Eq.  of  Motion  by  Displacement -based  Central  Difference  Method  * 

*  * 
************************************************************************ 

c 

subroutine  esolv  (time,nel,nnd,ndof ,  iaccnumout,  iforce,xmass, 
+  force, pint, rifix,d, v,a,xc,yc, mode, e, rkout, 

+  maxstp, dn, dp, inital,meq, sigmaP, sigmaN,epslonP, 

+  epslonN, PLalphaP, PLalphaN, PLrP, PLrN,elplas, 

+  elpout,proutv, si, s2, s3, s4, iplane) 

c 

implicit  real*8  (a-h,o-z) 

dimension  xmass (1) , force (1) ,pint (1) , mode (8,1) , rkout (3,1) 

dimension  d(l),v(l),a(l),xc(l),yc(l),e(12,l)  ,rifix(l) 

dimension  dn(l) , dp(l) ,ag (3) , elplas (1) , elpout (1) 

dimension  sigmaP(l) , sigmaN(l) ,epslonP(l) , epslonN (1) ,proutv(l) 

dimension  PLalphaP(l) ,PLalphaN(l) ,PLrP(l)  ,PLrN(l) 

dimension  si (1) , s2 (1) , s3 (1) , s4 (1) 

common  /box  1/  iprob, delta, alpha, toler, gravity 
common  /box  4/  npnts,numif ,nnaf ,ndisi,nveli, 
+  kfpnts(6) , jnode(30) , jdir(30)  ,  jaf (30) , 

+  ndnod(lO) ,kdis(10) ,nvnod(10) ,kvel(10) 

common  /box  4a/g, disi (10) , veli (10) , f (10)  ,  omega (10)  , 
+  ff (500,6) ,ta(500,6) , tf (500,6)  ,aa(500,6) 

c   define  initial  condition  that  time  is  within  the  range 
c     —  see  subroutine  finter  in  detail 
c 

noyes  =  0 
c 

do  100  i=l,meq 

d(i)  =  0.0 

v(i)  =  0.0 

a(i)  =  0.0 

100  forced)  =0.0 

c 

if  (inital  .eq.  0)  go  to  50 
if  (inital  .eq.  2)  go  to  60 
c 

c   set  index  no.  of  d,v,  and  put  initial  value  into  them 
c 

do  65  i=l,ndisi 

ind  =  (ndnod(i) -1) *ndof+kdis (i) 
d(ind)  =  disi (i) 
65  continue 
c 

if  (inital  .ne.  3)  go  to  50 
c 

60  continue 

do  70  i=l,nveli 

inv  =  (nvnod(i) -1) *ndof+kvel (i) 
v(inv)  =  veli(i) 
70  continue 
c 

c   compute  displacement  (nstep=0)  before  first  step  displ . (nstep=l) 
c  by  Central  Difference  Method 
c 

50  continue 
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do  75  i=l,meq 

dn(i)  =  d(i)-delta*v(i)+0.5*delta*delta*a(i) 
75  continue 
c 

nstep  =  0 

nskip  =  0 
c 
c 

160  do  120  i=l,meq 
120  forced)  =  0. 
c 

cc  impose  gravity  load 
c 

do  140  i=2,meq,2 
140  forced)  =  -1 .  *gravity*xmass  (i) 
c 

c@d — 20janl996 — cshih 
c 

cc  impose  impact  force 
c 

c   set  index  no.  of  nodal  d.o.f.  applied  by  impact  force 
c  compute  external  force  at  that  time  step  by  linear  interpolation 
c 

if  (iforce  .eq.  0)  go  to  167 
c 

do  165  n=l,nnaf 

imn  =  ( jnode (n) -1) *ndof+jdir (n) 

if  (iforce  .ne.  1)  go  to  163 

ma  =  jaf(n) 

mb  =  kfpnts(ma) 

call  finter  (f f (l,ma) , tf (l,ma) ,mb, time,pht,noyes) 

force (imn)  =  pht  +  force (imn) 

go  to  165 
c 

163  if (iforce.eq.2) force (imn) =f (n) *ds in (omega (n) *time) + force (imn) 

if (iforce . eq. 3 ) force (imn) =f (n) *dcos (omega (n) *time) +force (imn) 
165  continue 
c 

167  continue 
c 

cc   impose  ground  acceleration 
c 

if  (iacc  .eq.  0)  go  to  25 
c 

do  5  i=l,3 
5  agd)  =  0.0 
c 

c  compute  acceleration  at  that  time  step  by  linear  interpolation 
c 

if  (iacc  .It.  4)  go  to  13 

if  (iacc  .eq.  4)  go  to  12 
c 

do  11  j=l,3 

call  finter  (aa (1, j ) , ta (1, j ) ,npnts, time, pht, noyes) 

ag ( J )  =  Pht 

11  continue 
go  to  14 

c 

12  continue 

call  finter  (aa(l, 1) , ta (1, 1) ,npnts, time,pht, noyes) 
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ag(l)  =  pht 

call  f inter  (aa (1, 2) , ta (1, 2) , npnts, time, pht, noyes) 
ag(3)  =  pht 
go  to  14 
c 

13  call  f inter  (aa (1, 1) , ta (1, 1) ,npnts, time, pht, noyes) 
ag(iacc)  =  pht 

c 

c   compute  external  force  by  using  d'Alembert  principle  to  account  for 

c   inertia  force 

c 

14  do  20  i=l,nnd 
ii  =  (i-l)*ndof 
do  21  j=ii+l,ii+2 

force(j)  =  force (j ) -xmass (j ) *ag(j-ii) *g 
21  continue 
20  continue 
25  continue 
c 

c  initialized  index  for  average  stress  output 
c 
c 

do  998  ig  =  l,numout 
proutv(ig)  =  O.dO 
998  continue 
c 

do  999  io  =  l,nnd 
sl(io)  =  O.OdO 


c 


s2(io)  =  0. 

.  OdO 

s3(io)  =  0, 

,  OdO 

s4(io)  =  0, 

,  OdO 

999  continue 

do  991  ki  =  l,4*nel 

elplas(ki)  =  0 . OdO 

elpout(ki)  =  O.OdO 
991  continue 
c 
c   compute  element  internal  nodal  forces  (  pint  ) 


do  22  i=l,meq 

22  pint(i)=0. 

do  23  i=l,nel 

neo  =  int  (mode  (8,  i)  ) 

if  (neo. It. 0  .or.  neo.eq.ll)  goto  23 

call  fintiso8  (i,ndof , xc , yc , mode ,  e,d,pint,  sigmaP, sigmaN, 

+  epslonP,epslonN, PLalphaP, PLalphaN, PLrP, PLrN, 

+  elplas , si , s2 , s3 , s4 , iplane , time ) 

23  continue 

do  360  i=l,nel 
do  360  k=l,4 
n  =  (i-l)*4 
do  370  j=l,6 

m  =  (n+k-1) *6+j 

sigmaP(m)  =  sigmaN (m) 

epslonP(m)  =  epslonN(m) 

PLalphaP (m)  =  PLalphaN (m) 
370  continue 


76 


PLrP(n+k)  =  PLrN(n+k) 
360  continue 
c 

c   compute  displacement  of  next  time  step ~dp(j) 

c      and  velocity  &  acceleration  of  this  time  step  v(j)  &  a(j) 

c  by  Displacement -based  Central  Difference  Method 
c 

do  170  i=l,nnd 
m  =  (i-1) *ndof 
do  171  j=m+l,m+2 

if  (int(rifix(j) )  .eq.  1)  go  to  171 
if  (xmass(j)  .le.  0.0)  go  to  171 
al  =  2.0+alpha*delta 
a2  =  2.0-alpha*delta 
cl  =  (2.0/al)*delta*delta 
c2  =  4.0/al 
c3  =  a2/al 

c4  =  (force (j ) -pint (j )) /xmass (j ) 
dp(j)  =  cl*c4+c2*d(j)-c3*dn(j) 
v(j)  =  (dp(j)-dn(j))/(2.0*delta) 
a(j)  =  (dp(j)-2.0*d(j)+dn(j) )/(delta*delta) 
171  continue 
170  continue 
c 
c 

c  record  plastic  element  number 
c 

if  (nskip  .eq.  iprob)  then 
write (6,*)  'nstep=' ,nstep 
ij  =  0 

do  556  i=l,4*nel 

if  (elplas(i)  .It.  1.0)  go  to  556 
ij  =  ij  +  1 
elpout(ij)  =  elplas(i) 
556  continue 

write  (6,2021) 
2021  f ormat ( / , '   Plastic  element  no  [element  no. Gauss  point  no]  =') 
write  (6,2031)  (elpout (i) , i=l, i j ) 

2031  format (8 (3x,f 5.1)) 

if  (ij  .eq.  0)  write (6, 2032) 

2032  formate      NONE') 
endif 

c 

c   print  out  response  results  requested  every  "iprob"  step 

c 

if  (nskip  .ne.  iprob)  go  to  195 
if  (numout  .ne.  0)  then 

call  prout  (numout,  rkout,d, v,a, si, s2, s3 , s4,ndof ,nstep, 
+  time,proutv) 

nskip  =  0 
endif 
195  if  (nstep  .ge.  maxstp)  go  to  225 
c 

c   check  whether  the  increment  of  each  nodal  displacement  is  less  than 
c   the  tolerance  limit,  i.e.  TOLER. 
c 

if(nstep  .It.  200)  go  to  198 
ddmax  =  0 . OdO 
do  197  i=l,meq 

yn  =  dabs (dp(i) -d(i) ) 
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197  if  (ddmax  .It.  yn)  ddmax  =  yn 
if(ddmax  .It.  toler)  then 

write  (*,1100)  ddmax 
1100  format  (/,'  STOP  --  displ . increment  <  tolerance  limit ' , lpd9 .3 , /) 
go  to  225 
endif 

198  continue 
c 

do  200  i=l,meq 

dn(i)  =  d(i) 

d(i)   =  dp(i) 
200  continue 
c 

nskip  =  nskip  +  1 

nstep  =  nstep  +  1 

time   =  time  +  delta 

go  to  160 
c 

225  continue 
c 

c  record  plastic  element  numberfor  last  time  step 
c 

write (6,*)  ' nstep= ', nstep 

ij  =  0 

do  555  i=l,4*nel 

if  (elplas(i)  .It.  1.0)  go  to  555 

ij  =  ij  +  1 

elpout(ij)  =  elplas(i) 
555  continue 

write  (6,2001) 
2001  format(/,'   Plastic  element  no  =>[Element  no. Gauss  point  no]  =') 

write  (6,2011)  (elpout (i) , i=l, ij ) 
2011  format(8(3x,f5.1) ) 

if  (ij  .eq.  0)  write (6, 2022) 
2022  format ('      NONE') 

print  *, 'complete! ' 
c 

return 

end 
c 


************************************************************************ 

*  * 

*  calculate  each  time  step's  external  force  and  acceleration  by  linear* 

*  interpolation  because  of  mismatch  between  the  time  interval  of  * 

*  force  &  ace.  input  data,  and  the  time  increment  * 

*  * 

************************************************************************ 

c 

subroutine  f inter  (pp, w,kc, t,pht,noyes) 
c 

implicit  real*8  (a-h,o-z) 
dimension  pp(kc),w(kc) 
c 

if  (t.lt.w(l)  .or.  t.gt.w(kc))  then 
if  (noyes  .eq.  1)  go  to   900 
write (6, *) •  ' 
write (6,*) '  WARNING  --  time  out  of  the  range  of  time-force  pairs' 
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write (6,*)'  ZERO  FORCE  is  given  at  that  time.' 

write (6,*)'  the  out -of -the -range  time  is  '  ,t 

write (6, *) '  ' 
print  * ,  '  ' 

print  *,  '  WARNING  --  time  out  of  the  range  of  time-force  pairs' 
print  *,  '  ZERO  FORCE  is  given  at  that  time. ' 

print  *,  ■  the  out -of -the- range  time  is  ',t 

print  * ,  '  ' 
noyes  =  1 
900     pht  =0.0 
go  to  140 
endif 
c 

do  100  1=2, kc 
if  (  t  .le.  w(l))  go  to  120 
100  continue 
c 

120  pht  =  pp(l-l)  +  (t-w(l-l))*(pp(l)-pp(l-l))/(w(l)-w(l-l)) 
c 

140  return 
c 

c@d-22octl995-c . shih 
c 

end 

c 

••••••••••••••••••••••••••••a******************************************* 

*  * 

*  compute  element  internal  nodal  forces  (  pint  )  by  * 

*  4-node  plane  isoparametric  element  (  4-node  quadrilateral  element  )  * 

*  * 

•••••a****************************************************************** 

c 

subroutine  fintiso8  (i,ndof,xc,yc, mode, e,d, pint, sigmaP, sigmaN, 
+  epslonP,epslonN,PLalphaP,PLalphaN, PLrP, PLrN, elplas, 

+  si, s2, s3, s4, iplane, time) 

c 

implicit  real*8  (a-h,o-z) 

dimension  xc ( 1 ) , yc ( 1 ) , mode (8,l),e(12,l),d(l), pint ( 1 ) 

dimension  disl (8) , f (8) ,bmx(3,8) ,te(3,3) 

dimension  epsx(4) , epsy (4) , epsxy (4) 

dimension  ed(6, 6) , sigma (3) , sig(4, 3) , sigma3d(6) 

dimension  Depslon ( 6 ) , PLalpha ( 6 ) , elplas ( 1 ) 

dimension  sigmaP(l) , sigmaN (1) ,epslonP(l) ,epslonN(l) 

dimension  PLalphaP(l) ,PLalphaN(l) ,PLrP(l) ,PLrN(l) 

dimension  si ( 1 ) , s2 ( 1 ) , s3 ( 1 ) , s4 ( 1 ) 

common  /box  1/  iprob, delta, alpha, toler, gravity 

common  /Gauss/  s(4),t(4),w(4) 

data  phi/3.1415926535898/ 
c 

c   set  node  no .  for  element 
c 

nl  =  int  (mode  (1,  i) ) 

n2  =  int  (mode  (2,  i)  ) 

n3  =  int  ( mode  ( 3  ,  i )  ) 

n4  =  int  (mode  (4,  i)  ) 
c 

c  node (8)  =  1  element  condition  (  read  in  sub"readata"  ) 

c   if  node (8)  =  neo  =  11   >  element  has  a  negative  Jacobian 

c 

neo=  int  (mode  (8,  i)  ) 
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c 

c 
c 

c 
c 
c 


c 

c 
c 


if  (neo.eq.ll)  go  to  13 
set  index  no.  of  nodes 


nln  = 
n2n  = 

n3n  = 
n4n  = 


(nl-l)*ndof 
(n2-l)*ndof 
(n3-l)*ndof 
(n4-l)*ndof 


set  node  coordinates 


xl 
x2 
x3 

x4 

yi 
y2 
y3 

y4 


xc (nl) 
xc(n2) 
xc (n3 ) 
xc (n4) 
yc(nl) 
yc(n2) 
yc (n3 ) 
yc (n4) 


c 
c 
c 


set  node  displacement 


dlx 
dly 
d2x 
d2y 
d3x 
d3y 
d4x 
d4y 


d(nln+l) 
d(nln+2) 
d(n2n+l) 
d(n2n+2) 
d(n3n+l) 
d(n3n+2) 
d(n4n+l) 
d(n4n+2) 


defomed  coordinates 


xld  =  xl 
yld  =  yl 
x2d  =  x2 
y2d  =  y2 
x3d  =  x3 
y3d  =  y3 
x4d  =  x4 
y4d  =  y4 


dlx 
dly 
d2x 
d2y 
d3x 
d3y 
d4x 
d4y 


c 
c 
c 


check  distored  element 

A123  =  dabs(xld*y2d+x2d*y3d+x3d*yld-xld*y3d-x2d*yld-x3d*y2d) 
A134  =  dabs(xld*y3d+x3d*y4d+x4d*yld-xld*y4d-x3d*yld-x4d*y3d) 
A234  =  dabs (X2d*y3d+x3d*y4d+x4d*y2d-x2d*y4d-x3d*y2d-x4d*y3d) 
A124  =  dabs (Xld*y2d+x2d*y4d+x4d*yld-xld*y4d-x2d*yld-x4d*y2d) 

if  (  A123  .It.  1.0d-07  )  go  to  8 

if  (  A134  .It.  1.0d-07  )  go  to  8 

if  (  A234  .It.  1.0d-07  )  go  to  8 

if  (  A124  .It.  1.0d-07  )  go  to  8 

if  (  dabs(A123+A134-A234-A124)  .gt.  1.0d-07  )  go  to  8 

go  to  5 


8  write (6, 10)  i 
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10  format (lx, ' --  Distorted  element  no. 
mode  (8,  i)  =11.0 
go  to  13 


M5) 


c 
c 
c 
c 
c 


5  continue 

find  rotational  angle  of  element, 

based  on  the  lines  composite  with  centroid  and  each  nodal  point 

the  coordinates  of  the  centroid  of  element  =  (  xcd,  ycd  ) 

xcdu  =  (xl  +  x2  +  x3  +  x4)  /  4.d0 

ycdu  =  (yl  +  y2  +  y3  +  y4)  /  4.d0 

xcdd  =  (xld  +  x2d  +  x3d  +  x4d)  /  4.d0 

ycdd  =  (yld  +  y2d  +  y3d  +  y4d)  /  4.d0 


50 


60 


call  thetafind  (xcdu/ycdu/xcdd/ycdd,xld,yld,xl,yl, thetal) 

call  thetafind  (xcdu, ycdu, xcdd,ycdd, x2d, y2d, x2, y2, theta2) 

call  thetafind  (xcdu, ycdu, xcdd, ycdd, x3d, y3d, x3, y3, theta3) 

call  thetafind  (xcdu, ycdu, xcdd, ycdd, x4d, y4d, x4, y4, theta4) 

theta  =  (thetal  +  theta2  +  theta3  +  theta4)  /  4.d0 

if  (abs (theta)  .le.  1.0d-5)  theta  =  0 . OdO 

if  (theta  .It.  2.0d+0*phi)  go  to  60 

theta  =  theta  -  2.0d+0*phi 

go  to  50 

continue 


c 
c 
c 


c 
c 

) 

c 


form  transformation  matrix  global  to  local 


ell  = 
cl2  = 
c21  = 
c22  = 


dcos  ( 
dsin( 
-cl2 
ell 


theta 
theta 


compute  local  node  coord,  and  displ.  (  rotation  from  global  to  local 


xll  = 
yll  = 
xl2  = 
yl2  = 
xl3  = 
yl3  = 
xl4  = 
yl4  = 


cll*xl 
c21*xl 
cll*x2 
c21*x2 
cll*x3 
c21*x3 
cll*x4 
c21*x4 


cl2*yl 
c22*yl 
cl2*y2 
c22*y2 
cl2*y3 
c22*y3 
Cl2*y4 
c22*y4 


disl(l) 
disl(2) 
disl(3) 
disl(4) 
disl(5) 
disl(6) 
disl(7) 
disl(8) 

do  100  j=l, 
epsx(j) 
epsy(j) 
epsxy(j)  = 
f(2*j-l)  = 
100  f(2*j   )  = 


cll*dlx 
c21*dlx 
cll*d2x 
c21*d2x 
cll*d3x 
c21*d3x 
cll*d4x 
c21*d4x 


cl2*dly 
c22*dly 
cl2*d2y 
c22*d2y 
cl2*d3y 
c22*d3y 
cl2*d4y 
c22*d4y 


OdO 
OdO 
OdO 
O.OdO 
O.OdO 
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c 
c 
c 

c 
c 


c 
c 
c 


c 
c 
c 


(1 

dO-t(k) ) 

/ 

4 

.dO 

(1 

dO-t(k)) 

/ 

4 

dO 

(1 

dO+t(k)) 

/ 

4 

dO 

(1 

dO+t(k) ) 

/ 

4 

dO 

coeff.  of  4  point  Gauss  Quadrature  have  been  transmitted  from 

common  /Gauss/  s (4) , t (4) ,w{4)    in  sub."bmass" 

set  shape  function  Ni 


do  200  k=l,4 
shpfl  =  (l.dO-s(k) ) 
shpf2  =  (l.d0+s(k) ) 
shpf3  =  (l.d0+s(k) ) 
shpf4  =  (l.dO-s(k) ) 


derivatives  of  shape  function  w.r.t.  s,t 

slds  =  -(l.dO-t(k) )  /  4.d0 

s2ds  =   (l.dO-t(k) )  /  4.d0 

s3ds  =   (l.d0+t(k) )  /  4.d0 

s4ds  =  -(l.d0+t(k))  /  4. dO 

sldt  =  -(l.dO-s(k))  /  4. dO 

s2dt  =  -(l.d0+s(k))  /  4.d0 

s3dt  =   (l.d0+s(k) )  /  4.d0 

s4dt  =   (l.dO-s(k) )  /  4. dO 

derivatives  of  global  cood.  x,y  w.r.t.  s,t 

xs  =  slds*xll  +  s2ds*xl2  +  s3ds*xl3  +  s4ds*xl4 
ys  =  slds*yll  +  s2ds*yl2  +  s3ds*yl3  +  s4ds*yl4 
xt  =  sldt*xll  +  s2dt*xl2  +  s3dt*xl3  +  s4dt*xl4 
yt  =  sldt*yll  +  s2dt*yl2  +  s3dt*yl3  +  s4dt*yl4 

detJ  =  xs*yt  -  ys*xt 


if  (detJ  .le.  O.dO)  then 


c 
c 
c 


print  * , 
write (6, * ) 
write (6, *) 
write (6, *) 
write (6, *) 
mode  ( 8 ,  i ) 
go  to  13 
endif 

form  B  matrix 


Negative  Jacobian,  ele.  no  =  ' ,  i 

Negative  Jacobian,  ele.  no  =  '  ,  i 

node  no .  1,2,3,4  =  ' ,  nl , n2 , n3 , n4 
local  x=  ' ,xll,xl2,xl3,xl4, 'global  x= 
local  y=  ' ,yll,yl2,yl3,yl4, 'global  y= 
=  11.0 


J= 
J= 


detJ 
detJ 


,  xl ,  x2  ,  x3  ,  x4 
,  yl ,  y2 ,  y3 ,  y4 


blx 
bly 
b2x 
b2y 
b3x 
b3y 
b4x 
b4y 


yt*slds 
xs*sldt 
yt*s2ds 
xs*s2dt 
yt*s3ds 
xs*s3dt 
yt*s4ds 
xs*s4dt 


ys*sldt 
xt*slds 
ys*s2dt 
xt*s2ds 
ys*s3dt 
xt*s3ds 
ys*s4dt 
xt*s4ds 


bmx(l,l) 
bmx (1,2) 
bmx (1,3) 
bmx (1,4) 
bmx (1,5) 
bmx (1,6) 
bmx (1,7) 
bmx (1,8) 


blx 

O.dO 

b2x 

O.dO 

b3x 

O.dO 

b4x 

O.dO 


82 


bmx(2,l) 
bmx (2, 2) 
bmx (2,3) 
bmx(2,4) 
bmx (2,5) 
bmx (2,6) 
bmx (2,7) 
bmx (2, 8) 
bmx(3,l) 
bmx (3, 2) 
bmx (3, 3) 
bmx (3, 4) 
bmx (3, 5) 
bmx (3,6) 
bmx (3,7) 
bmx(3,8) 


O.dO 

bly 

O.dO 

b2y 

O.dO 

b3y 

O.dO 

b4y 

bly 

blx 

b2y 

b2x 

b3y 

b3x 

b4y 

b4x 


c   compute  local  element  strain 

c 

do  300  j=l,8 
epsx(k)   =  epsx(k)   + 
epsy(k)   =  epsy(k)   + 
epsxy(k)  =  epsxy(k)  + 
300  continue 


bmx(l, j)*disl(j)  /  detJ 
bmx (2, j)*disl(j)  /  detJ 
bmx (3, j)*disl(j)  /  detJ 


assign  value  of  new  epslon  (epslonN(..) 
check  plane  strain  problem 


if  (iplane  .eq.  2)  then 
m  =  (  (i-l)*4+k-l)*6 
epslonN(m+l) 
epslonN (m+2) 
epslonN(m+3) 


epslonN (m+4) 
epslonN (m+5) 
epslonN (m+6) 
end  if 


epsx(k) 

epsy(k) 

O.OdO 

epsxy (k) 

O.OdO 

O.OdO 


c 
c 
c 


check  plane  stress  problem 


if  (iplane  .eq.  1)  then 
m  =  (  (i-l)*4+k-l)*6 
epslonN (m+1) 


epslonN (m+2) 
epslonN (m+3) 
epslonN (m+4) 
epslonN (m+5) 
epslonN (m+6) 
end  if 


epsx(k) 

epsy(k) 

- (poisson/ (1-poisson) ) * (epsx(k) +epsy (k) 

epsxy (k) 

O.OdO 

O.OdO 


c 
c 

c 


c 
c 

c 


calculate  delta  epslon 

do  352  in  =  1,6 
Depslon(in)  =  epslonN (m+in)  -  epslonP(m+in) 
352  continue 

material  properties  for  element  (  b=thickness,  eyng=  E  ) 
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mtyp    =  int  ( mode  ( 5 ,  i )  ) 

mtyp2   =  int (e(l,intyp)  ) 

eyng    =  e(3,mtyp) 

if  (eyng  .le.  O.OdO)  goto  13 

poisson  =  e (4, mtyp) 

ft      =  e(5,mtyp) 

coh     =  e ( 6 , mtyp ) 

phiangle  =  e (7, mtyp) *phi/180 .dO 

kickPL    =  int (e(8,mtyp) ) 

eyngt    =  e(9,mtyp) 

beta     =  e(10,mtyp) 

b        =  e( 11, mtyp) 

tau      =  e( 12, mtyp) 
c 
c   call  constitutive  coefficient  from  subroutine. 


c 


call  elastd3d  (eyng, poisson, ed) 


if  (kickPL  .eq.  1)  then 

do  402  m  =  1, 6 

do  502  n  =  1,6 

ed(m,n)  =  exp(-time/tau) *ed(m,n) 
502  continue 
402  continue 

end  if 
c 

c  calculate  trial  stress (sigma  =  stresses  at  gauss  pts) 
c 

ny  =  (i-1) *4  +  k 
PLr  =  PLrP(ny) 
do  353  m=l,6 
ny  =  ( (i-l)*4+k-l)*6+m 
sigma3d(m)  =  sigmaP(ny) 
PLalpha(m)  =  PLalphaP(ny) 
353  continue 

do  400  m  =  1,6 
do  500  n  =  1,6 
sigma3d(m)  =  sigma3d(m)  +  ed(m,n) *Depslon(n) 
500  continue 
400  continue 
c 

sigmean  =  sigma3d(l) +sigma3d(2) +sigma3d(3) /3 . OdO 
c 

c  apply  plasticity  model  :  Drucker-Prager  yield  criterion 
c 

if (kickPL  .ne.  1) 
+call  PLmodel (sigma3d, eyng, poisson, ft, coh, phiangle, kickPL, 
+  eyngt , beta , PLalpha , PLr , sigmean , mtyp2 , elplas , i , k) 

c 

PLrN( (i-l)*4+k)  =  PLr 
do  800  m=l, 6 
ny  =  ( (i-l)*4+k-l)*6+m 
sigmaN(ny)  =  sigma3d(m) 
PLalphaN(ny)  =  PLalpha (m) 
800  continue 
c 

c   transfer  stresses  from  3d  to  plane  strain  2d  condition 
c 

sigma(l)  =  sigma3d(l) 
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c 
c 
c 
c 


sigma (2)  =  sigma3d(2) 
sigma(3)  =  sigma3d(4) 

transform  local  nodal  stresses  into  global 

--sig(i,k)  =  sigma{i=x- (1) ,y- (2) ,xy- (3) ;  gauss  pts(k=l,4) 

te(l,l)  =  cll*cll 

ted, 2)  =  cl2*cl2 

te(l,3)  =  cll*cl2 

te(2,l)  =  te(l, 2) 

te(2,2)  =  te(l,l) 

te(2,3)  =  -ted, 3) 

te(3,l)  =  -2.0d0*te(l,3) 

te(3,2)  =   2.0d0*te(l,3) 

te(3,3)  =  te(l,l)  -  ted, 2) 


810 
777 


do  777  mi=l,3 
sig(k,mi)  =  0. 
do  810  ni=l,3 
sig(k,mi)  = 
continue 
continue 


OdO 


sig(k,mi)  +  te (ni, mi )* sigma (ni) 


c 
c 
c 


compute  local  element  internal  nodal  forces 


do  600  m=l,8 
do  700  n=l,3 
f(m)  =  f(m)  + 

7  00  continue 
600  continue 


b*w(k) *bmx(n,m) * sigma (n) 


c 
c 
c 
c 

c 

c 
c 


200  continue 

Extrapolate  stress  from  Gauss's  pts  to  node 
and  calculate  average  nodal  stress 

call  stress (sig,nl,n2,n3,n4, si, s2, s3, s4) 

transform  local  nodal  forces  to  global 


c 
c 
c 


fix  = 

cll*f  (1 

+ 

c21*f (2) 

fly  = 

cl2*f (1 

+ 

c22*f  (2) 

f2x  = 

cll*f (3' 

+ 

c21*f (4) 

f2y  = 

cl2*f (3' 

+ 

c22*f  (4) 

f3x  = 

cll*f (5 

+ 

c21*f (6) 

f3y  = 

cl2*f (5 

+ 

c22*f  (6) 

f4x  = 

cll*f (7 

+ 

c21*f (8) 

f4y  = 

cl2*f (7 

+ 

c22*f (8) 

assemble  to  global  nodal  force  matrix 


pint  (internal) 


pint 
pint 
pint 
pint 
pint 
pint 
pint 
pint 


(nln+1) 
(nln+2) 
(n2n+l) 
(n2n+2) 
(n3n+l) 
(n3n+2) 
(n4n+l) 
(n4n+2) 


pint 
pint 
pint 
pint 
pint 
pint 
pint 
pint 


(nln+1 
(nln+2 
(n2n+l 
(n2n+2 
(n3n+l 
(n3n+2 
(n4n+l 
(n4n+2 


fix 
fly 
f2x 
f2y 
f3x 
f3y 
f4x 
f4y 
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c 

13  return 
end 

c 

************************************************************************ 

*  Drucker-Prager  yield  criterion  with  non-associate  flow  rule  and      * 

*  linear  combination  of  isotropic  and  kinematic  hardening/softening    * 

*  Krieq  and  Key's  radial-return  algorithm  for  elastoplastic  case       * 

*  (  kickPL  =  3  for  Drucker-Prager)  * 

*  (  beta  =  0.  for  kinematic  &  =1  for  isotropic  hardening)      * 
************************************************************************ 

c 

subroutine  PLmodel (sigma3d, eyng,poisson, ft, coh, phiangle, kickPL, 
+  eyngt,beta, PLalpha, PLr, sigmean/mtyp2 , elplas,  i,k) 

c 

implicit  real*8  (a-h,o-z) 

dimension  PLalpha (1) , sigma3d(l) ,xi (6) , elplas (1) 
c 

Ep  =  eyngt/(1.0d0  -  eyngt/eyng) 
gshear  =  eyng/ (2 . OdO* (1 . OdO  +  poisson) ) 
cl  =  2.0d0*gshear  +  Ep*2 . OdO/3 . OdO 
c2  =  2.0dO/3.0dO*beta*Ep 
c3  =  2.0dO/3.0dO*(1.0dO-beta)*Ep 
c4  =  2.0d0*gshear 
c 

c  calculate  xi-trial 
c 

do  10  j=l,6 
xi(j)  =  sigma3d(j)  -  PLalpha(j) 
10  continue 
c 

c  calculate  deviatoric  part  of  xi 
c 

ximean  =  (xi (l)+xi (2) +xi (3) ) /3 . OdO 
c 

c  calculate  the  radius  of  cross  section  of  cone/cylinder 
c 

if  ( (mtyp2  .eq.  2 ) .and. (kickPL  .eq.  2))  phiangle  =  0 . OdO 
if(kickPL  .eq.  2)  then 
b  =  O.OdO 

ak  =  ft/dsqrt(3.0d0) 
endif 

if  (kickPL  .eq.  3)  then 
b  =  dsqrt (12 .d0 ) *dsin (phiangle) / (3 . OdO+dsin (phiangle) ) 
ak  =  dsqrt (12 .d0) * (-coh) *dcos (phiangle) / (3 .dO+dsin (phiangle) ) 
endif 

r  =  dabs (dsqrt (2 .d0) * (ak+b*sigmean) ) 
if  (r  .gt.  PLr)  PLr  =  r 
c 

do  20  j  =1,3 
xi(j)  =  xi(j)  -  ximean 
20  continue 
c 

c  check  if  elastic 
c 

xxi  =  xi(l)*xi(l)+xi(2)*xi(2)+xi(3)*xi(3)+ 
+      2.d0*(xi(4)*xi(4)+xi(5)*xi(5)+xi(6)*xi(6) ) 
yn  =  PLr* PLr 

if  (xxi  .le.  yn)  go  to  100 
c 
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c  Plastic  phase:  calculate  unit  normal — N  (also  store  in  "xi") 
c 

xxi  =  dsqrt(xxi) 

do  30  j  =1,6 
xi  (j )  =  xi ( j) /xxi 
30  continue 
c 

c  calculate  lambda-tilde  ("A") 
c 

A  =  (xxi-PLr) /cl 
c 

c  update  stresses 
c 

PLr  =  PLr  +  c2*A 

tl  =  c3*A 

t2  =  c4*A 

do  40  j=l,6 
PLalpha(j)  =  PLalpha(j)  +  tl*xi(j) 
sigma3d(j)  =  sigma3d(j)  -  t2*xi(j) 
40  continue 
c 

c   record  the  element  number  if  plastic 
c 

j  =  4*(i-l)+k 

elplas(j)  =  dble(i)  +  0 .ld0*dble (k) 


100  return 

end 

c 
************************************************************************ 

*  * 

*  compute  element  internal  nodal  streses (sigma)  by  * 

*  8-node  solid  isoparametric  element  (8-node  hexahedral  element)       * 

*  * 

************************************************************************ 

c 

subroutine  stress (sig,nl,n2,n3,n4, si, s2,  s3,s4) 
c 

implicit  real*8  (a-h,o-z) 

dimension  x(4) ,y(4) 

dimension  sig (4,3) , signodeX (4) , signodeY(4)  ,  signodeXY(4) 

dimension  sN(4,4),sl(l),s2(l),s3(l),s4(l) 
c 

common  /box  1/  iprob, delta, alpha, toler, gravity 
c 

c  Interpolate  stresses  from  gauss  points  to  nodes 
c  (  si(idof,io)  ;  idof=dof,  io=  local  node  no.  ) 
c 

constant  =  dsqrt(3.0d0) 

x(l)  =  -constant 

y(l)  =  -constant 
c 

x(2)  =   constant 

y(2)  =  -constant 
c 

x(3)  =   constant 

y(3)  =  constant 
c 

x(4)  =  -constant 
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203 


y(4)  =   constant 

do  203  j  =  1,4 

sN(j,l)  =  (I.d0-x(j))*(l.d0-y(j))/4.d0 

sN(j,2)  =  (I.d0+x(j))*(l.d0-y(j))/4.d0 

sN(j,3)  =  (l.dO+x(j) )*(l.dO+y(j) )/4.d0 

sN(j,4)  =  (l.dO-x(j) )*(l.d0+y(j) )/4.d0 

continue 


do  501  j  =  1,4 

signodeX(j)  =  0 . OdO 
signodeY(j)  =  0 . OdO 
signodeXY(j)  =  0 . OdO 
do  502  k  =  1,4 

signodeX(j)  =  signodeX(j)  ■* 
signodeY(j)  =  signodeY(j)  + 
signodeXY(j)  =  signodeXY(j) 
502  continue 
501  continue 

sl(nl)  =  sl(nl)  +  signodeX(l) 

sl(n2)  =  sl(n2)  +  signodeX(2) 

sl(n3)  =  sl(n3)  +  signodeX(3) 

sl(n4)  =  sl(n4)  +  signodeX(4) 


sN(j,k)*sig(k,l) 
sN(j,k)*sig(k,2) 
+  sN(j,k)*sig(k,3) 


s2(nl)  =  s2(nl)  +  signodeY(l) 

s2(n2)  =  s2(n2)  +  signodeY(2) 

s2(n3)  =  s2(n3)  +  signodeY(3) 

s2(n4)  =  s2(n4)  +  signodeY(4) 


s3(nl)  =  s3(nl)  +  signodeXY(l) 

s3(n2)  =  s3(n2)  +  signodeXY(2) 

s3(n3)  =  s3(n3)  +  signodeXY(3) 

s3(n4)  =  s3(n4)  +  signodeXY(4) 
c 

s4(nl)  =  s4(nl)  +  1 . OdO 

s4(n2)  =  s4(n2)  +  1 . OdO 

s4(n3)  =  s4(n3)  +  1 . OdO 

s4(n4)  =  s4(n4)  +  1 . OdO 
13  return 

end 
************************************************************************ 

*  * 

*  print  out  final  results  requested  at  specific  points  &  directions    * 

*  * 
************************************************************************ 


subroutine  prout  (numout, rkout,d, v,a, si, s2, s3 , s4,ndof , 
+  nstep, time,proutv) 

c 

c  This  subroutine  controls  output  of  displacement,  velocity, 
c  acceleration,  and  stresses, 
c 

implicit  real*8  (a-h,o-z) 

dimension  d(ndof , 1) , v(ndof , 1) ,a(ndof , 1) 

dimension  rkout(3,l) ,proutv(l) ,sl (1) , s2 (1) , s3 (1)  ,  s4 (1) 
c 

do  100  i2=l, numout 
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node  =  int (rkout (1, i2) ) 
idva  =  int (rkout (2 , i2) ) 
idof  =  int (rkout (3, i2) ) 

if  (idva  .eg.  0)  proutv(i2)  =  d( idof, node) 

if  (idva  .eq.  1)  proutv(i2)  =  v(idof,node) 

if  (idva  .eq.  2)  proutv(i2)  =  a (idof, node) 

if  (idva  .eq.  3)  then 

if  (idof  .eq.  1)  proutv(i2)  =  si (node) /s4 (node) 
if  (idof  .eq.  2)  proutv(i2)  =  s2 (node) /s4 (node) 
if  (idof  .eq.  3)  proutv(i2)  =  s3 (node) /s4 (node) 
end  if 
100  continue 

write  (*,200)  nstep 
200  format ( '  [ '  ,i7,  '  ]  ' ) 

write (7, 300)  time, (proutv(i) , i=l,numout) 
300  format (lx, 'time  =' ,el2 . 5, 6 (lx,e9 .3) /19x, 6 (lx,e9.3) / 

+  19x,6(lx,e9.3)/19x,6(lx,e9.3)/19x,6(lx,e9.3)/19x,6(lx,e9.3) / 
+  19x,6(lx,e9.3)/19x,6(lx,e9.3)/19x,6(lx,e9.3)/19x,6(lx,e9.3)/ 
+  19x,6(lx,e9.3)/19x,6(lx,e9.3)/19x,6(lx,e9.3)/19x,6(lx,e9.3) / 
+  19x,6(lx,e9.3)/19x,6(lx,e9.3)/19x,6(lx,e9.3)/19x,6(lx,e9.3) / 
+  19x,6(lx,e9.3)/19x,6(lx,e9.3)/19x,6(lx,e9.3)/19x,6(lx,e9.3) / 
+  19x,6(lx,e9.3)/19x,6(lx,e9.3)/19x,6(lx,e9.3)/19x,6(lx,e9.3)/ 
+  19x,6(lx,e9.3)/19x,6(lx,e9.3)/19x,6(lx,e9.3)/19x,6(lx,e9.3)/ 
+  19x,6(lx,e9.3)/19x,6(lx,e9.3)/19x,6(lx,e9.3)/19x,6(lx,e9.3)/ 
+  19x,6(lx,e9.3)/19x,6(lx,e9.3)/19x,6(lx,e9.3)/19x,6(lx,e9.3)/ 
+  19x,6(lx,e9.3)/19x,6(lx,e9.3)/19x,6(lx,e9.3)/19x,6(lx,e9.3)/ 
+  19x,6(lx,e9.3)/19x,6(lx,e9.3)/19x,6(lx,e9.3)/19x,6(lx,e9.3)/ 
+  19x,6(lx,e9.3)/19x,6(lx,e9.3)/19x,6(lx,e9.3)/19x,6(lx,e9.3)/ 
+  19x,6(lx,e9.3) /19x, 6 (lx, e9 . 3 ) /19x, 6 (lx, e9 . 3) /19x, 6 (lx,e9 .3) / 
+       19x,6(lx,e9.3) ) 

return 
end 


************************************************************************ 

*  * 

*  read  in  node  &  element  data  ,  material  properties  * 

*  ace.  &  force  data,  initial  condition  data,  and  output  request   * 

*  * 

•••••A****************************************************************** 

c 

subroutine  readata  (nnd,nel,nummat,numout, iacc,ndof , mode, 
+  iforce, imesh,xc,yc, rif ix, e, rkout, inital) 

c 

implicit  real*8  (a-h,o-z) 

dimension  rif ix(l) ,rnode (8,1) , rkout (3,1) ,xc (1) ,yc (1) , e (12 , 1) 

dimension  if ixx(2) ,node(8) , kout (3) 

c 

common  /box  1/  iprob, delta, alpha, toler, gravity- 
common  /box  4/  npnts,numif ,nnaf ,ndisi,nveli , 
+  kfpnts(6) , jnode(30) , jdir(30) , jaf (30) , 
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+  ndnod(lO) ,kdis(10) ,nvnod(10)  ,kvel(10) 

common  /box  4a/g,disi (10) , veli (10) , f (10) , omega (10) , 
+  ff (500, 6), ta (500, 6) , tf (500,6) ,aa(500,6) 

c 

meq  =  nnd*ndof 
c 

do  100  i=l,meq 
100  rifix(i)  =  0.0 
c 

c   read  &  write  nodal  coord.  &   B.C. 
c   skipped  nodes  will  be  automatically  generated 
c 

write(6,550) 

1  =  0 
110  read(5,*)   n,xc (n) ,yc (n) , (ifixx(i) , i=l, 2) 
c 

do  120  j=l,2 

nl  =  (n-l)*ndof 

rifix(nl+j)  =  dble (if ixx( j ) ) 
120  continue 
c 

nl  =  1+1 
c 

if  (1  .eq.  0)  go  to  130 
c 

dl  =  n-1 

dxl  =  (xc(n)-xc(l))/dl 

dyl  =  (yc(n)-yc(l))/dl 
c 

130  continue 

I  =  1+1 
c 

if  (n-1)  180,170,150 
c 

150  xc(l)  =  xc(l-l)+dxl 

yc(l)  =  yc(l-l)+dyl 
c 

do  155  j=l,2 

II  =  (l-l)*ndof 
12  =  (l-2)*ndof 
rifix(ll+j)  =  rifix(12+j) 

155  continue 
c 

if  (imesh  .eq.  1)  then 
write(6,570)  l,xc(l) ,yc(l) , (int (rifix(i) ) , i=ll+l, 11+2) 
endif 
go  to  130 
c 
170   continue 

if  (imesh  .eq.  1)  then 

write (6, 570)  n,xc (n) ,yc (n) , (int (rifix(k) ) ,k=nl+l,nl+2) 
endif 

if  (nnd-n)  180,190,110 
c 

180  continue 

write(6,580)  n 
stop 
c 

190  continue 
c 
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c  read  &  write  element  connectivity 
c 

write(6,600) 

1  =  0 
200  read(5, *)   n,node(l) ,node(2) ,node(3) ,node(4) ,node(5) ,node(6) , 
+  node ( 7 )  , node ( 8 ) , knt 

c 

do  212  i=l,8 

rnode(i,n)  =  dble(node (i) ) 
212  continue 
c 

nl  =  1+1 

if  (knt  .eq.  0)  knt  =  1 
245  continue 

1  =  1  +  1 

if  (n-1)  260,255,250 
c 
250   rnode(l,l)  =  rnode(l,l-l)  +  dble(knt) 

mode  (2,1)  =  mode  (2, 1-1)  +  dble(knt) 

mode  (3,1)  =  mode  (3, 1-1)  +  dble(knt) 

rnode(4,l)  =  rnode(4,l-l)  +  dble(knt) 

rnode(5,l)  =  mode(5,l-l) 

mode  (6,1)  =  mode  (6,1-1) 

mode(7,l)  =  rnode(7,l-l)  +  dble(knt) 

mode  (8,1)  =  mode  (8, 1-1) 

go  to  245 
c 

255  continue 

if  (imesh  .eq.  1)  then 

write  (6,  620)  (k,  int  (mode  (1,  k)  )  ,  int  (mode  (2  ,  k)  )  ,  int  (mode  (3  ,  k)  )  , 
+  int  (mode  (4,k)  )  ,  int  (mode  (5,k)  )  ,  int  (mode  (6,k)  )  , 

+  int  (mode  (7,k)  )  ,  int  (mode  (8,k)  )  ,  k=nl,n) 

endif 

if  (nel-n)  260,270,200 
c 

2  60  continue 

write(6,630)  n 

stop 
c 

270  continue 
c 

c  read  &  write  material  properties 
c 

do  300  i=l,nummat 

read (5,*)   k,e (l,k) , e (2,k) ,e (3,k) ,e (4,k) , e (5,k) 

read (5,*)   e(6, k) ,e (7, k) ,e(8, k) ,e (9,k) , e (10, k) , e (11, k) ,e (12 ,k) 

write(6,960) 

write (6, 970)  k,int(e(l,k) ) , e (2 ,k) , e (3 ,k) , e (4,k) , e (5,  k) 

write(6,980) 

write (6, 990)  e (6,k) ,e(7,k) ,int(e(8,k) ) , e (9,k) ,e (10, k)  ,  e (11, k) 
300  continue 
c 

c  read  &  write  no.  of  acceleration  history,  and  data  of  time-acc.  pairs 
c   in  each  ace.  history 
c 

if  (iacc  .eq.  0)  go  to  760 
c 

read ( 5 , * )     nacc , npnts 

read ( 5 , * )     g 

write (6, 1170) 
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write (6, 1180)  nacc,npnts 
c     write(6,1185) 
c     write(6,1186)  g 
c 

do  721  i=l,nacc 

read (5,*)   ( ta ( j , i ) , aa  ( j , i )  ,  j  =1 , npnts ) 
721  continue 
c 

760  continue 

if  (iforce  .eq.  0)  go  to  762 

if  (iforce  .ne.  1)  go  to  763 
c 

c   read  &  write  no.  of  impact  force  history,  no.  of  nodes,  and  data  of 
c   time-acc.  pairs  in  each  ace.  history 


c 


c 
c 


read (5,*)     numif,nnaf 

write (6, 1300) 

write (6, 1310)numif ,nnaf 

if  (numif  .eq.  0)  go  to  76! 


do  755  i=l, numif 

read (5,*)   kfpnts(i) 

jf  =  kfpnts (i) 

read (5,*)   (tf ( j , i) , f f ( j , i) , j=l, jf ) 

write(6,1295) 

do  785  j=l,jf 

write (6, 1210)  i, j , tf ( j , i) , f f ( j , i) 
785  continue 
755  continue 
c 

768  continue 
write(6,1315) 

c 

c   read  &  write  data  of  position  applied  by  arbitrary  shape  impact  force 

c   function  (  node,  d.o.f.,  history  no.  ) 

c 

do  769  i=l,nnaf 

read (5, *)     jnode (i) , jdir (i) , jaf (i) 

write (6, 122 5) jnode(i) , jdir (i)  ,  jaf (i) 

769  continue 
go  to  762 

c 

c   read  &  write  data  of  position  applied  by  sinusoidal  impact  force 

c   function  (  node,  d.o.f.,  history  no.  ) 

c 

763  continue 

write (6, 1316) 

read (5,*)   nnaf 
c 

do  764  i=l,nnaf 

read ( 5 ,  * )   j  node ( i ) , j  dir ( i ) , f ( i ) , omega ( i ) 

write (6, 122  6)  jnode (i) , jdir (i) , f (i) , omega (i) 
7  64  continue 
c 

7  62  continue 

if  (inital  .eq.  0)  go  to  765 

if  (inital  .eq.  2)  go  to  782 
c 
c  read  &  write  data  of  displacement  type  I.C.  (  node,  dof,  I. C. value  ) 
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c 

read (5, *)   ndisi 

write(6,784) 

do  786  i=l, ndisi 

read (5, *)     ndnod(i) ,kdis (i)  ,disi (i) 

write (6, 783)  ndnod(i) ,kdis(i)  ,disi(i) 

786  continue 
c 

if  (inital  .ne.  3)  go  to  765 
c 

c   read  &  write  data  of  velocity  type  I.e.  (  node,  dof,  I. C. value  ) 
c 

782  read (5,*)   nveli 
write(6,787) 

do  788  i=l, nveli 

read (5,*)     nvnod(i) ,kvel (i)  , veli (i) 
write (6, 783)  nvnod(i) ,kvel(i) ,veli (i) 
788  continue 
c 

c   read  &  write  data  of  output  record  requested  (node,  d-v-a-type,  dof) 
c 

765  continue 

if  (numout  .eq.  0)  go  to  780 
c 

write(6,1228) 
write(7,1228) 
c 

do  77  0  i=l, numout 
read(5,*)  (kout ( j ) , j=l,  3) 
write(6,1229)  i, (kout ( j ) , j=l, 3) 
write(7,1229)  i, (kout ( j ) , j=l, 3) 
do  747  j=l,3 

rkout(j.i)  =  dble(kout(j) ) 
747  continue 
770  continue 
c 

7  80  continue 
c 
c 

550  format  (/2x, 'card  4',7x, 'nodal  point  data',/, 

+        lOx,  'node  no. ' ,3x, 'x-ordinate' ,2x,  'y-ordinate' , 
+        4x, 'ifx' ,4x, 'ify' ) 
570  format  (lOx, i5 , 2x, 2f 12 . 3 , 2i8) 
580  format  (17x, 'nodal  point  error  n  =',i5) 
c 

600  format (/2x, 'card   5     element  data ' / , lOx, ' ele.  no. ' , lx, 'node-1 ' , 
+         lx,  'node-2 ' , lx,  'node-3 ' , lx,  'node-4 ' , lx,  'mat-typ' , 
+         '  row-no' , '  col-no' , '  ele-cond. ' ) 
620  format  (lOx, i6, lx, i6, lx, i6, lx, i6, lx, i6, lx, i6, lx, i6, lx, i6, lx, i6) 
630  format  (17x, ' element  number  error  n  =',i5) 
c 

783  format  (i8, 6x, i6, lOx, fl2 . 5) 

784  format  (5x, 'node ', 5x, ' disp-dir' , 8x, ' initial  disp',/) 

787  format  (5x, 'node' , 5x, 'vel-dir' , 8x, ' initial  vel',/) 
c 

1170  format  (/,2x, 'card   12 ', 6x, 'prescribed  acceleration  card'/) 

1180  format  (lOx, 'total  no.  of  acceleration  histories  =', 

+    i5, /, lOx, ' total  no.  of  time-acc.  pairs  in  each  ace.  history  =', 

+    15) 
c 
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960  format  (/2x,'card   6  &  7    material  property  data' , /, 10x, 

+ 'material   material   mass      Youngs     Poisson    tensile ' /10x, 
+' group  no.  type  no.   density   modulus    ratio    strength') 


970  format 

980  format 

+ ' angle 

990  format 

c 

1200  format 

+ 
1210  format 
1220  format 

+ 
1225  format 
122  6  format 

1228  format 
+  'seq. 

1229  format 
c@d 

1295  format 

+ 
c 
13  00  format 
1310  format 

+ 

1315  format 
+ 

1316  format 
+ 

c 

return 
end 


Ilx,i3,7x,i3,3x/el0.4,lx,el0.4,2x,f5.3,3x,el0.4) 

19x, 'cohesion   phi     yield     tangent  hardening' , /3 Ox 

criterion  modulus    rule     thick(b)') 
17x,el0.4,2x,f6.2,5x,i2,5x,e9.4,2x,f5.3,4x,f5.3) 

/,2x, 'card  13 ', 6x, 'acceleration-history  card',/, 

17x, ' pair  no . ' , 8x, ' time ' , 9x, ' ace ' / ) 

20x,i6,llx, i6,4x,el2 . 4,4x,el2 .4) 

/,2x, 'card  14 ', 6x, 'nodal  accele.  information  card',/, 

17x, '  node','   dir','   ace'/) 

15x,i5,8x,i5,14x,i5) 

15x,i5,8x,i5,14x,el0.4,3x,f5.3) 

/,2x, 'card  21   stress  output  information  card',/,14x, 

node*    d-(0),v-(l),a-(2),sig-(3)     x(l) ,y (2) ,xy (3 ) ' ) 
12x, i4, 6x, i4, 13x, i4, 22x, i4) 

/2x, 'card  12  &  13     impact  force  history  card',/,17x, 
force  history  no. '5x, 'pair  no. ' , 5x, ' time' , 9x, ' iforce ' ) 

/2x, 'card  11  prescribed  impact  force') 

20x, 'total  no.  of  impact  force  history          =',i5,/, 

20x, 'total  no.  of  nodes  applied  by  impact  force  =',i5) 

/2x, 'card  14  nodal  impact  force  information',/, 

15x, 'node  no.  x-(l),y(2)    force  history  no.') 

/2x, 'card   14  sinusoidal  force  information',/, 

15x, 'node  no.  x- (1) ,y- (2) , z- (3)     ampli.      freg. ' ) 


************************************************************************ 

*  * 

*  use  arc-cos  function  to  find  theta  value  * 

*  theta  =  dacos (  numerator  /  denominator  )  * 

*  * 

************************************************************************ 

c 

subroutine  thetafind  (xcdl,ycdl, xcd2 ,ycd2 ,xd, yd, xu,yu, theta) 
c 

implicit  real*8  (a-h,o-z) 

data  phi/3.1415926535898/ 
c 

xau  =  xu  -  xcdl 

yau  =  yu  -  ycdl 


xad  =  xd  -  xcd2 

yad  =  yd  -  ycd2 
c 
c 

c  compute  the  length  of  undeformed  shape 
c 

rLu  =  dsqrt (xau*xau  +  yau*yau) 
c 

c  compute  the  length  of  deformed  shape 
c 

rLd  =  dsqrt (xad*xad  +  yad*yad) 
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c  compute  the  unit  components  of  undeformed  vectors 
c 

euxl  =  xau/rLu 

euyl  =  yau/rLu 
c 

c  compute  the  unit  components  of  deformed  vectors 
c 

edxl  =  xad/rLd 

edyl  =  yad/rLd 
c 

c  compute  the  rigid  body  rotation 
c  xy  plane  rotation (thetaZ) 
c 

Cz  =  euxl* edyl  -  edxl* euyl 

theta  =  dasin(Cz) 

if  ( (xu*xd  .It.  O.dO) .and. (yu*yd  .It.  O.dO))  then 
if  (thetaZ  .ge.  O.dO)  thetaZ  =  phi  -  thetaZ 
if  (thetaZ  .It.  O.dO)  thetaZ  =  -phi  -  thetaZ 

endif 

return 

end 
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Solid3D  Source  Code 
(S3DP.FOR) 
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************************************************************************ 

*  * 

*  Program  for  Dynamic  of  3  Dimensional  Solid  Element  Problem         * 

*  * 

************************************************************************ 

*  * 

*  This  program  is  for  analysis  of  3-D  dynamic  problem  and  developed   * 

*  on  the  basis  of  :  * 

*  (1)  Traditionally  Co-Rotational  Approach  * 

*  (2)  Explicit  Time  Integration  Method  (central  difference)  * 

*  (3)  Lumped  Mass  Modeling  * 

*  (4)  8-node  Solid  Isoparametric  Element  * 

*  (5)  8-point  Gaussian  Integration  * 

*  (6)  Elastic-lenear  work-hardening  Model  * 

*  (7)  Viscoelastic  Model  * 

*  (8)  von  Mises  yield  criterion  * 

*  (9)  Drucker-Prager  yield  criterion  * 

*  * 

*  Tatsana  Nilaward  * 

*  Purdue  University  * 

*  04/10/96 

*  * 
************************************************************************ 

c 

program  solid3D 

implicit  real*8  (a-h,o-z) 

dimension  ar (200000) 

maxq  =  200000 
c 

c   read  data  file  and  create  the  input  and  output  filenames 
c  (  extension  part  of  filename :_.dat,  .in,  and  .out  ,will  be  created 
automatically  ) 
c 

open  (5,  file= ' s3dp.dat ' ) 

open  (6,  file= ' s3dp.in' ) 

open  (7,  file= ' s3dp.out ' ) 
c 

call  allocate  (ar,maxq) 
c 

print  *,  'TOTALLY  COMPLETE  !  ' 

close  (unit=5) 

close  (unit=6) 

close  (unit=7) 
c 

stop 

end 

c 

************************************************************************ 

*  * 

*  set  index  no.  of  each  variable  by  allocation  method  * 

*  read  in  control  parameters   &  call  main  subroutines  * 

*  * 

************************************************************************ 

c 

subroutine  allocate  (ar,maxq) 
c 

implicit  real*8  (a-h,o-z) 

dimension  ar(maxq) 

dimension  head (20) 

common  /box  1/  iprob, delta, alpha, toler, gravity 
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common  /box  5/  maxcyc ,  maxout ,  maxmat 
c 

do  100  i=l,maxq 
100  ar(i)  =  O.dO 
c 

read (5, 130)   head 

write (6 ,140)  head 

write(6,150) 
c 

c   iprob  is  number  of  time-steps  skipped  between  output 

c 

read ( 5 ,  * )     iprob, nnd, nel , nummat , numout , ndof , maxstp , delta , alpha 

*  ,  toler,  gravity- 
write  (6,160)  iprob, nnd, nel, nummat, numout, ndof, maxstp, delta, alpha 

*  ,  toler,  gravity- 
read  (5, *)     iacc, i force, inital, imesh 
write (6, 170)  iacc, i force, inital, imesh 

c 

c  allocation  of  ar()  

c 

meq    =  nnd  *  ndof 

maxel   =  nel 

maxmat  =  nummat 

maxnod  =  nnd 

npint   =  1 

nxmass  =  npint  +  meq 

na     =  nxmass  +  meq 

nv     =  na  +  meq 

nd     =  nv  +  meq 

nxc     =  nd  +  meq 

nyc    =  nxc  +  maxnod 

nzc    =  nyc  +  maxnod 

nforce  =  nzc  +  maxnod 

nifix  =  nforce  +  meq 
c 

c   element  node  connectivity 
c 

nnode   =  nifix+meq 
c 

c  material  properties 
c 

nem    =  nnode+12*maxel 
c 
c   output  record 


c 


maxcyc  =  maxstp 
maxout  =  numout 

if  (maxstp  .eq.  0)  maxcyc  =  1 
if  (numout  .eq.  0)  maxout  =  1 

nkout   =  nem  +  ll*maxmat 
ndp    =  nkout  +  3*maxout 
ndn    =  ndp  +  meq 

nsigmaP  =  ndn  +  meq 
nsigmaN  =  nsigmaP  +  8*6*nel 
nepslonP  =  nsigmaN  +  8*6*nel 
nepslonN  =  nepslonP  +  8*6*nel 
nPLalphaP  =  nepslonN  +  8*6*nel 
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nPLalphaN  =  nPLalphaP  +  8*6*nel 
nPLrP  =  nPLalphaN  +  8*6*nel 
nPLrN  =  nPLrP  +  8*nel 
nelplas  =  nPLrN  +8*nel 
nelpout  =  nelplas  +8*nel 
nsl  =  nelpout  +  8*nel 


ns2 

= 

nsl 

+ 

nnd 

ns3 

= 

ns2 

+ 

nnd 

ns4 

= 

ns3 

+ 

nnd 

ns5 

= 

ns4 

+ 

nnd 

ns6 

= 

ns5 

+ 

nnd 

ns7 

= 

ns6 

+ 

nnd 

nproutv  =  ns7  +  nnd 

maxindex  =  nproutv  +  maxout 
c 

if  (  maxindex  .gt.  maxq  )  then 

print  *,  '  There  is  not  enough  dimension  available.  ' 

print  * ,  '  Please  increase  no .  in  ar ( . . . )  &  maxq= ' , maxindex 

stop 

endif 
c 

c  call  subroutines 
c 

call  readata  (nnd, nel, nummat, numout, iacc, ndof , ar (nnode) , iforce, 
+  imesh,ar (nxc) , ar (nyc) , ar (nzc) ,ar (nifix) , ar (nem) , 

+  ar (nkout) , inital) 

print  *,  'call  readata  --  complete!' 
c 

call  bmass  (nel , nnd, ndof , ar (nem) ,ar (nnode) ,ar (nxc) ,ar (nyc) , 
+  ar (nzc) ,ar (nxmass) ) 

print  *,  'call  bmass  complete!' 


time 


O.dO 


call  esolv  (time, nel, nnd, ndof , iacc,numout, iforce, ar (nxmass) , 

+  ar(nforce) ,ar(npint) , ar (nifix) ,ar(nd) ,ar(nv) ,ar(na) , 

+  ar (nxc)  ,  ar (nyc) , ar (nzc) , ar (nnode) , ar (nem) , ar (nkout) , 

+  maxstp, ar (ndn) , ar (ndp) , inital , meq, ar (nsigmaP) , 

+  ar (nsigmaN) , ar (nepslonP) , ar (nepslonN) , ar (nPLalphaP) , 

+  ar (nPLalphaN) ,ar(nPLrP) ,ar(nPLrN) ,ar(nelplas) , 

+  ar(nelpout) ,ar(nproutv) ,ar(nsl) ,ar(ns2) ,ar(ns3) ,ar(ns4) 

+  ar (ns5) , ar (ns6 ) , ar (ns7) ) 


print 


'call  esolv  complete! 


c 

c 

130  format 

(20a4) 

140  format 

(/2x, 'card   l',5x,20a4) 

150  format 

(lx,80('-')) 

160  format 

(2x,  'card  2 ', 5x,  'parameter 

card ' , / , 

+ 

15x, "no  of  time-steps  skipped  between 

outputs 

+ 

15x, 'number  of  nodes 

= 

,il0,/, 

+ 

15x, 'number  of  elements 

= 

,il0,/, 

+ 

15x, 'number  of  materials 

= 

,il0,/, 

+ 

15x, 'number  of  output  req 

= 

,il0,/, 

+ 

15x, 'no.  of  d.o.f/node 

= 

,il0,/, 

+ 

15x, 'no.  of  time  steps 

= 

,il0,/, 

+ 

15x, 'time  increment 

= 

,el0.3, 

/, 

+ 

15x, 'coeff  of  mass  damping 

= 

,el0.3. 

/, 

+ 

15x, ' tolerance  limit 

= 

,el0.3, 

/, 

=',i6,/, 
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+  15x, 'acceleration  of  gravity  =',fl0.5,/) 

170  format  (2x, 'card  3 ' , 5x, ' index  card',/, 

+  15x, ' index  for  accel .       =',il0,/, 

+  15x, 'index  for  force       =',il0,/, 

+  15x, 'index  for  I.  C.       =',il0,/, 

+  15x, 'index  for  mesh  output (1)  or  not(0)  =',i4) 

c 

return 

end 
************************************************************************ 

*  * 

*  calaulate  translational  mass  of  all  nodes  in  x,  y,  and  z-dir        * 

*  and  form  Mass-matrix  all  at  one  time  * 

*  * 
************************************************************************ 

c 

subroutine  bmass  (nel,nnd,ndof ,e,rnode,xc,yc, zc,xmass) 
c 

implicit  real*8  (a-h,o-z) 

dimension  mode (12, 1) ,xc(l) ,yc(l) ,zc(l) ,e(ll,l) ,xmass (1) 
c 

c   calculate  mass  for  each  node 
c 

meq  =  nnd*ndof 
c 

do  10  i=l,meq 
10  xmass(i)  =  O.dO 
c 

c   loop  through  all  element 
c 

do  20  i=l,nel 
c 

c   set  node  no.  and  node  coord,  for  element 
c 

nl  =  int  (mode  (1,  i) 

n2  =  int  ( mode  ( 2 ,  i ) 

n3  =  int  (mode  (3  ,  i) 

n4  =  int  (mode  (4,  i) 

n5  =  int  (mode  (5,  i) 

n6  =  int  (mode  (6,  i) 

n7  =  int  (mode  (7,  i) 

n8  =  int  (mode  (8,  i) 


xl  =  xc (nl) 
x2  =  xc(n2) 
x3  =  xc(n3) 
x4  =  xc(n4) 
x5  =  xc (n5) 
x6  =  xc(n6) 
x7  =  xc(n7) 
x8  =  xc(n8) 

yl  =  yc(nl) 
y2  =  yc(n2) 
y3  =  yc (n3) 
y4  =  yc (n4) 
y5  =  yc(n5) 
y6  =  yc (n6) 
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y7  =  yc(n7) 

y8  =  yc (n8) 
c 

zl  =  zc(nl) 

z2  =  zc(n2) 

z3  =  zc(n3) 

z4  =  zc(n4) 

z5  =  zc(n5) 

z6  =  zc(n6) 

z7  =  zc(n7) 

z8  =  zc(n8) 
c 

c  material  properties  for  element  (  b=thickness,  rho=mass  density  ) 
c 

mtyp  =  int  (mode  (9,  i) ) 

rho  =  e ( 2 , mtyp ) 
c 

c   calculate  area  of  base  triangle  and  height 
c   aijk  =  area  of  triangle  with  nodes  i,j,and  k  are  vertex 
c 

al25  =  dabs(0.5d0*( (x2*z5+x5*zl+xl*z2) - (x2*zl+x5*z2+xl*z5) ) ) 

a438  =  dabs(0.5d0* ( (x4*z8+x8*z3+x3*z4) - (x4*z3+x8*z4+x3*z8) ) ) 

a265  =  dabs(0.5d0* ( (x5*z6+x6*z2+x2*z5) - (x5*z2+x6*z5+x2*z6) ) ) 

a378  =  dabs(0.5d0*( (x7*z8+x8*z3+x3*z7) - (x7*z3+x8*z7+x3*z8) ) ) 

hi  =  dabs(y4-yl) 

h2  =  dabs(y3-y2) 

h3  =  dabs(y7-y6) 

h4  =  dabs(y8-y5) 
c 

c   calculate  average  area  and  height 
c 

aavgl  =  (al25+a438) 12  .d0 

havgl  =  (hl+h2+h4)/3.d0 

aavg2  =  (a265+a378) 12  .dO 

havg2  =  (h2+h3+h4) /3 -dO 


c 

c 

compute 

each  half  mass  of  element  i  th 

c 

rmhalfl  =  rho* aavgl* havgl 
rmhalf2  =  rho*aavg2*havg2 
totmass  =  rmhalfl +rmhalf 2 

c 

c 

calculate  index  of  mass  storage 

c 

ml  = 

(nl-l)*ndof 

m2  = 

(n2-l)*ndof 

m3  = 

(n3-l)*ndof 

m4  = 

(n4-l)*ndof 

m5  = 

(n5-l)*ndof 

m6  = 

(n6-l)*ndof 

m7  = 

(n7-l)*ndof 

c 
c 

m8  = 

(n8-l)*ndof 

assemble 

i  to  global  mass  matrix  xmass 

(8) 

c 

xmass 

(ml+1)  =  xmass (ml+1)  +  totmass/8 

.dO 

xmass 

(ml+2)  =  xmass (ml+2)  +  totmass/8 

.dO 

xmass 

(ml+3)  =  xmass (ml+3)  +  totmass/8 

.dO 

xmass 

(m2+l)  =  xmass (m2+l)  +  totmass/8 

.dO 

xmass 

;(m2+2)  =  xmass (m2+2)  +  totmass/8 

.dO 
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xmass (m2+3 
xmass (m3+l 
xmass (m3+2 
xmass (m3+3 
xmass (m4+l 
xmass (m4+2 
xmass (m4+3 
xmass (m5+l 
xmass (m5+2 
xmass (m5+3 
xmass (m6+l 
xmass (m6+2 
xmass (m6+3 
xmass (m7+l 
xmass (m7+2 
xmass (m7+3 
xmass (m8+l 
xmass (m8+2 
xmass (m8+3 
continue 
return 
end 


xmass 
xmass 
xmass 
xmass 
xmass 
xmass 
xmass 
xmass 
xmass 
xmass 
xmass 
xmass 
xmass 
xmass 
xmass 
xmass 
xmass 
xmass 
xmass 


(m2+3 
(m3+l 
(m3+2 
(m3+3 
(m4+l 
(m4+2 
(m4+3 
(m5+l 
(m5+2 
(m5+3 
(m6+l 
(m6+2 
(m6+3 
(m7+l 
(m7+2 
(m7+3 
(m8+l 
(m8+2 
(m8+3 


totmass 
totmass 
totmass 
totmass 
totmass 
totmass 
totmass 
totmass 
totmass 
totmass 
totmass 
totmass 
totmass 
totmass 
totmass 
totmass 
totmass 
totmass 
totmass 


/8.d0 
/8.d0 
/8.d0 
/8.d0 
/8.d0 
/8.d0 
/8.d0 
/8.d0 
/8.d0 
/8.d0 
/8.d0 
/8.d0 
/8.d0 
/8.d0 
/8.d0 
/8.d0 
/8.d0 
/8.d0 
/8.d0 


************************************************************************ 


*   compute  elastic 

* 


constitutive  coefficients  under   3D  condition 


************************************************************************ 


c 
c 


subroutine  elastd3d  (eyng,poisson, ed) 

implicit  real*8  (a-h,o-z) 
dimens ion  ed ( 6 , 6 ) 


poisson*eyng/ ( (1 .dO+poisson) * (l.dO-2 .d0*poisson) ) 


+  dla 


do  50  i=l 

,6 

do  60  j=l 

,6 

ed(i,j) 

=  O.dO 

60  continue 

50  continue 

G  =  eyng/ (2.d0*( 

dla   =  poisson*e}< 

ed(l,l)  = 

2.dO*G 

ed(2,2)  = 

ed(l,l) 

ed(3,3)  = 

ed ( 1 , 1 ) 

ed(4,4)  = 

G 

ed(5,5)  = 

ed(4,4) 

ed(6,6)  = 

ed(4,4) 

ed(l,2)  = 

dla 

ed(2,l)  = 

ed(l,2) 

ed(3,l)  = 

ed(l,2) 

ed(3,2)  = 

ed ( 1 , 2 ) 

ed(l,3)  = 

ed(l,2) 

ed(2,3)  = 

ed(l,2) 

return 

end 
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************************************************************************ 

*  * 

*  Solve  Eq.  of  Motion  by  Displacement -based  Central  Difference  Method  * 

*  * 
••••••••••A************************************************************* 

c 

subroutine  esolv  (time,nel,nnd,ndof , iacc,numout, iforce,xmass, 
+  force, pint, rif ix,d, v,a,xc,yc, zc,rnode,e,rkout, 

+  maxstp,dn,dp, inital,meq, sigmaP, sigmaN,epslonP, 

+  epslonN,PLalphaP,PLalphaN,PLrP, PLrN,elplas, 

+  elpout , proutv, si , s2 , s3 , s4 , s5 , s6 , s7 ) 

c 

implicit  real*8  (a-h,o-z) 

dimension  xmass (1) , force (1) ,pint (1) , mode (12 , 1) , rkout (3,1) 
dimension  d(l),v(l),a(l),xc(l),yc(l),zc(l),e(ll,l) ,rifix(l) 
dimension  dn ( 1 ) , dp ( 1 ) , ag ( 3 ) , elplas ( 1 ) , elpout ( 1 ) 
dimension  sigmaP(l) ,sigmaN(l) ,epslonP(l) ,epslonN(l) , proutv (1) 
dimension  PLalphaP ( 1 ) , PLalphaN ( 1 ) , PLrP ( 1 ) , PLrN ( 1 ) 
dimension  si  (1)  ,  s2  (1)  ,  s3  (1)  ,  s4  (1)  ,  s5  (1)  ,  s6  (1)  ,  s7  (1) 
c 

common  /box  1/  iprob, delta, alpha, toler, gravity 
common  /box  4/  npnts,numif ,nnaf ,ndisi,nveli, 
+  kfpnts(6) , jnode(30) , jdir(30) , jaf (30) , 

+  ndnod(lO) ,kdis(10) ,nvnod(10) ,kvel (10) 

common  /box  4a/g,disi (10) , veli (10) , f (10) , omega (10) , 
+  ff (500,6) ,ta(500,6) ,tf (500,6) ,aa(500,6) 

c 

c   define  initial  condition  that  time  is  within  the  range 
c     --  see  subroutine  finter  in  detail 

noyes  =  0 
c 

do  100  i=l,meq 

d(i)  =  O.dO 

v(i)  =  O.dO 

a(i)  =  O.dO 

100  continue 

c 

if  (inital  .eq.  0)  go  to  50 
if  (inital  .eq.  2)  go  to  60 
c 

c   set  index  no.  of  d,v,  and  put  initial  value  into  them 
c 

do  65  i=l,ndisi 

ind  =  (ndnod(i) -1) *ndof+kdis (i) 
d(ind)  =  disi (i) 
65  continue 

if  (inital  .ne.  3)  go  to  50 

60  continue 

do  70  i=l,nveli 

inv  =  (nvnod(i) -1) *ndof+kvel (i) 

v(inv)  =  veli (i) 
70  continue 

c 

c   compute  displacement  (nstep=0)  before  first  step  displ . (nstep=l) 

c   by  Central  Difference  Method 

c 

50  continue 


c 
c 
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do  75  i=l,meq 
dn(i)  =  d(i)-delta*v(i)+0.5*delta*delta*a(i) 
75  continue 
c 

nstep  =  0 
nskip  =  0 
c 

160  continue 
c 

do  120  i=l,meq 
forced)  =  O.dO 
120  continue 
c 

c  impose  gravity  load 
c 

do  140  i=2,meq,3 
forced)  =  -1  .d0*gravity*xmass  (i) 
140  continue 
c 

c  impose  impact  force 

c   set  index  no.  of  nodal  d.o.f.  applied  by  impact  force 
c  compute  external  force  at  that  time  step  by  linear  interpolation 
c 

if  (iforce  .eq.  0)  go  to  167 

do  165  n=l,nnaf 
imn  =  ( jnode(n) -1) *ndof+jdir (n) 
if  (iforce  .ne.  1)  go  to  163 
ma  =  jaf(n) 
mb  =  kfpnts(ma) 

call  finter  (f f (l,ma) , tf (l,ma) ,mb, time, pht, noyes) 
force (imn)  =  pht  +  force (imn) 
go  to  165 
c 

163  if (iforce. eq. 2) force (imn) =f (n) *ds in (omega (n) *time) + force (imn) 
if (iforce.eq.3) force (imn) =f (n) *dcos (omega (n) *time) + force (imn) 
165  continue 
167  continue 

c 

c  impose  ground  acceleration 
c 

if  (iacc  .eq.  0)  go  to  25 
c 

do  5  i=l,3 
ag(i)  =  O.dO 
5  continue 
c 

c   compute  acceleration  at  that  time  step  by  linear  interpolation 
c 

if  (iacc  .It.  4)  go  to  13 
if  (iacc  .eq.  4)  go  to  12 
c 

do  11  j=l,3 
call  finter  (aa(l, j ) , ta (1, j ) , npnts, time, pht ,noyes) 
ag ( j )  =  pht 

11  continue 
go  to  14 

c 

12  continue 

call  finter  (aa (1, 1) , ta (1, 1) ,npnts, time, pht, noyes) 


104 


ag(l)  =  pht 

call  f inter  (aa (1, 2) , ta (1, 2) ,npnts, time, pht, noyes) 
ag ( 3 )  =  pht 
go  to  14 
c 

13  call  f inter  (aa(l, 1) , ta (1, 1) ,npnts, time, pht, noyes) 
ag(iacc)  =  pht 

c 

c  compute  external  force  by  using  d'Alembert  principle  to  account  for 

c   inertia  force 

c 

14  continue 

do  20  i=l,nnd 

ii  =  (i-l)*ndof 
do  21  j=ii+l,ii+3 

force(j)  =  force (j ) -xmass (j ) *ag(j-ii) *g 

21  continue 
20  continue 
2  5  continue 

c 

c  initialized  index  for  average  stress  output 

c 

c 

do  998  ig  =  l,numout 

proutv ( ig )  =  0 . dO 

998  continue 
c 

do  999  io  =  l,nnd 
sl(io)  =  O.OdO 
s2(io)  =  O.OdO 
s3(io)  =  O.OdO 
s4(io)  =  O.OdO 
s5(io)  =  O.OdO 
s6(io)  =  O.OdO 
s7(io)  =  O.OdO 

999  continue 
c 

do  991  ki  =  l,8*nel 
elplas(ki)  =  O.OdO 
elpout(ki)  =  O.OdO 
991  continue 
c 

c  compute  element  internal  nodal  forces (pint) 
c 

do  22  ix=l,meq 
pint(ix)  =  O.dO 

22  continue 
c 

do  23  i  =  l,nel 
neo  =  int  (mode  (12,  i)  ) 
if  (neo. It. 0   .or.  neo.eq.ll)  goto  23 
call  fintiso8  (i,ndof ,xc,yc, zc, mode, e, d, pint, sigmaP, 
+  sigmaN,epslonP, epslonN, PLalphaP, PLalphaN, 

+  PLrP, PLrN,elplas, si, s2, s3, s4, s5,s6, s7, time) 

23  continue 
c 

do  360  i=l,nel 

n  =  (i-l)*8 

do  360  k=l,8 
do  370  j=l,6 
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m  =  (n+k-1) *6+j 
sigmaP(m)  =  sigmaN(m) 
epslonP(m)  =  epslonN(m) 
PLalphaP(m)  =  PLalphaN(m) 
370  continue 

PLrP(n+k)  =  PLrN(n+k) 
360  continue 
c 

c  compute  displacement  of  next  time  step        dp(j) 

c  and  velocity  &  acceleration  of  this  time  step  v(j)  &  a(j) 

c  by  Displacement-based  Central  Difference  Method 
c 

do  170  i=l,nnd 

m  =  (i-1) *ndof 
do  171  j=m+l,m+3 
if  (int(rifix(j) )  .eq.  1)  go  to  171 
if  (xmass(j)  .le.  O.dO)  go  to  171 
al  =  2.d0+alpha*delta 
a2  =  2.d0-alpha*delta 
cl  =  (2.d0/al)*delta*delta 
c2  =  4.d0/al 
c3  =  a2/al 

c4  =  (force (j ) -pint (j )) /xmass (j ) 
dp(j)  =  Cl*c4+c2*d(j)-c3*dn(j) 
v(j)  =  (dp(j)-dn(j))/(2.d0*delta) 
a(j)  =  (dp(j)-2.d0*d(j)+dn(j) )/(delta*delta) 
171  continue 
170  continue 

c 

c 

c  record  plastic  element  number 

c 

if  (nskip  .eq.  iprob)  then 

write (6,*)  'nstep= ' ,nstep 

ij  =  0 

do  556  i=l,8*nel 

if  (elplas(i)  .It.  1.0)  go  to  556 

ij  =  ij  +  1 

elpout(ij)  =  elplas(i) 
556  continue 

write  (6,2021) 
2021  format(/, '   Plastic  element  no  [element  no. Gauss  point  no]  =') 

write  (6,2031)  (elpout (i) , i=l, ij ) 

2031  format(8(3x,f5.1) ) 

if  (ij  .eq.  0)  write (6, 2032) 

2032  format ('      NONE*) 
endif 

c 

c  print  out  response  results  requested  every  "iprob"  step 

c 

if  (nskip  .ne.  iprob)  go  to  195 

if  (numout  .ne.  0)  then 

call  prout  (numout, rkout,d, v,a, si, s2, s3, s4, s5, s6,  s7,ndof  ,nstep, 
+  time,proutv) 

nskip  =  0 

endif 
195  continue 

if  (nstep  .ge.  maxstp)  go  to  225 
c 
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c   check  whether  the  increment  of  each  nodal  displacement  is  less  than 

c   the  tolerance  limit,  i.e.  TOLER. 

c 

if(nstep  .It.  200)  go  to  198 

ddmax  =  0 . OdO 

do  197  i=l,meq 

yn  =  dabs (dp(i) -d(i) ) 

197  if  (ddmax  .It.  yn)  ddmax  =  yn 
if (ddmax  .It.  toler)  then 

write  (*,1100)  ddmax 
1100  format  (/,'  STOP  —  displ . increment  <  tolerance  limit ' , lpd9 .3 , /) 
go  to  225 
endif 

198  continue 
c 

do  200  i=l,meq 
dn(i)  =  d(i) 
d(i)   =  dp(i) 
200  continue 
c 

nskip  =  nskip  +  1 
nstep  =  nstep  +  1 
time   =  time  +  delta 
go  to  160 
c 

225  continue 
c 

c  record  plastic  element  numberfor  last  time  step 
c 

write (6,*)  'nstep= ', nstep 
ij  =  0 

do  555  i=l,8*nel 

if  (elplas(i)  .It.  1.0)  go  to  555 
ij  =  ij  +  1 

elpout(ij)  =  elplas(i) 
555  continue 

write  (6,2001) 
2001  format(/,'   Plastic  element  no  =>[Element  no. Gauss  point  no]  =  ') 

write  (6,2011)  (elpout (i) , i=l, ij ) 
2011  format(8(3x,f5.1) ) 

if  (ij  .eg.  0)  write (6, 2022) 
2022  format ( '      NONE' ) 
print  * , ' complete ! ' 
c 

return 
end 

c 

************************************************************************ 

*  * 

*  calculate  each  time  step's  external  force  and  acceleration  by  linear* 

*  interpolation  because  of  mismatch  between  the  time  interval  of  * 

*  force  &  ace.  input  data,  and  the  time  increment  * 

*  * 

************************************************************************ 
c 

subroutine  finter  (pp, w,kc, t,pht,noyes) 
c 

implicit  real*8  (a-h,o-z) 

dimension  pp(kc)  ,  w(kc) 
c 
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900 


if  (t.lt.w(l)  .or.  t.gt.w(kc))  then 

if  (noyes  .eq.  1)  go  to   900 

write(6, * 

write (6 

write (6 

write (6 

write (6 

print  * 

print  * 

print  * 

print  * 

print  * 

noyes  =  1 
pht  =0.0 

go  to  140 
end  if 


WARNING  —  time  out  of  the  range  of  time-force  pairs' 
ZERO  FORCE  is  given  at  that  time. ' 
the  out-of-the-range  time  is  ',t 


WARNING  —  time  out  of  the  range  of  time-force  pairs' 
ZERO  FORCE  is  given  at  that  time. ' 
the  out-of-the-range  time  is  ■ ,t 


do  100  1=2, kc 

if  (  t  .le.  w(l))  go  to  120 

100  continue 

c 

120  pht  =  pp(l-l)  +  (t-w(l-l))*(pp(l)-pp(l-l))/(w(l)-w(l-l)) 

c 

140  return 

c 

end 

c 
************************************************************************ 

*  * 

*  Compute  Element  Internal  Nodal  Forces  (pint)    by  * 

*  8-node  bicubic  isoparametric  element  (8-node  solid  element)         * 

*  * 

************************************************************************ 

c 

subroutine  fintiso8  (i,ndof ,xc,yc,  zc,  mode,  e,  d,  pint,  sigrnaP, 
+        sigmaN,epslonP,epslonN, PLalphaP, PLalphaN, PLrP, PLrN, elplas, 
+        si, s2, s3, s4, s5, s6, s7, time) 
c 

implicit  real*8  (a-h,o-z) 

dimension  xc(l)  ,yc(l)  ,  zc (1)  ,  mode (12, 1) ,e(ll, 1) ,d(l) ,pint (1) 

dimension  disl (24),f(24),b(6,24),te(6,6) 

dimension  epsx ( 8 ) , epsy ( 8 ) , epsz ( 8 ) , epsxy ( 8 ) , epsyz ( 8 ) , epszx ( 8 ) 

dimension  ed ( 6 , 6 ) , sigma ( 6 ) , sig (8,6) 

dimension  Depslon ( 6 )  , PLalpha ( 6 ) , elplas ( 1 ) 

dimension  sigmaP(l) ,sigmaN(l) ,epslonP(l) ,epslonN(l) 

dimension  PLalphaP ( 1 ) , PLalphaN ( 1 ) , PLrP ( 1 ) , PLrN ( 1 ) 

dimension  s(8),t(8),r(8),w(8) 

dimension  si  (1)  ,  s2  (1)  ,  s3  (1)  ,  s4  (1)  ,  s5  (1)  ,  s6  (1)  ,  s7  (1) 


c 
c 

c 


common  /box  1/  iprob, delta, alpha, toler, gravity 
data  phi/3. 141592653589793d0/ 

matching  local  to  global  nodal  number  for  one  element 


nl  =  int  (mode  (1,  i) 

n2  =  int  ( mode  ( 2 ,  i ) 

n3  =  int  ( mode  ( 3  ,  i ) 

n4  =  int  ( mode  ( 4 ,  i ) 

n5  =  int  (mode  (5,  i) 

n6  =  int  (mode  (6,  i) 


108 


n7  =  int  (mode  (7,  i)  ) 

n8  =  int  (mode  (8,  i)  ) 
c 

c  matching  local  to  global  coordinates  for  one  element 
c 

xl  =  xc(nl) 

x2  =  xc  (n2 ) 

x3  =  xc(n3) 

x4  =  xc(n4) 

x5  =  xc(n5) 

x6  =  xc(n6) 

x7  =  xc (n7) 

x8  =  xc(n8) 
c 

yl  =  yc(nl) 

y2  =  yc(n2) 

y3  =  yc (n3 ) 

y4  =  yc(n4) 

y5  =  yc(n5) 

y6  =  yc(n6) 

y7  =  yc(n7) 

y8  =  yc(n8) 
c 

zl  =  zc(nl) 

z2  =  zc(n2) 

z3  =  zc(n3) 

z4  =  zc(n4) 

z5  =  zc(n5) 

z6  =  zc(n6) 

z7  =  zc(n7) 

z8  =  zc(n8) 
c 

c  node (12)  =  1  element  condition  (  read  in  sub"readata"  ) 

c   if  node (12)  =  neo  =  11   >  element  has  a  negative  Jacobian 

c 

neo  =  int  (mode (12,  i)  ) 
c 

if  (neo.eq.ll)  go  to  13 
c 

c   set  index  no.  of  nodes 
c 

ml  =  (nl-l)*ndof 

m2  =  (n2-l)*ndof 

m3  =  (n3-l)*ndof 

m4  =  (n4-l)*ndof 

m5  =  (n5-l)*ndof 

m6  =  (n6-l)*ndof 

m7  =  (n7-l)*ndof 

m8  =  (n8-l)*ndof 
c 

c   set  node  displacement 
c 

dlx  =  d(ml+l) 

dly  =  d(ml+2) 

dlz  =  d(ml+3) 
c 

d2x  =  d(m2+l) 

d2y  =  d(m2+2) 

d2z  =  d(m2+3) 
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d3x  =  d(m3+l) 

d3y  =  d(m3+2) 

d3z  =  d(m3+3) 
c 

d4x  =  d(m4+l) 

d4y  =  d(m4+2) 

d4z  =  d(m4+3) 
c 

d5x  =  d(m5+l) 

d5y  =  d(m5+2) 

d5z  =  d(m5+3) 
c 

d6x  =  d(m6+l) 

d6y  =  d(m6+2) 

d6z  =  d(m6+3) 
c 

d7x  =  d(m7+l) 

d7y  =  d(m7+2) 

d7z  =  d(m7+3) 
c 

d8x  =  d(m8+l) 

d8y  =  d(m8+2) 

d8z  =  d(m8+3) 
c 

c  deformed  coordinates 
c 

xld  =  xl  +  dlx 

yld  =  yl  +  dly 

zld  =  zl  +  dlz 

x2d  =  x2  +  d2x 

y2d  =  y2  +  d2y 

z2d  =  z2  +  d2z 

x3d  =  x3  +  d3x 

y3d  =  y3  +  d3y 

z3d  =  z3  +  d3z 

x4d  =  x4  +  d4x 

y4d  =  y4  +  d4y 

z4d  =  z4  +  d4z 

x5d  =  x5  +  d5x 

y5d  =  y5  +  d5y 

z5d  =  z5  +  d5z 

x6d  =  x6  +  d6x 

y6d  =  y6  +  d6y 

z6d  =  z6  +  d6z 

x7d  =  x7  +  d7x 

y7d  =  y7  +  d7y 

z7d  =  z7  +  d7z 

x8d  =  x8  +  d8x 

y8d  =  y8  +  d8y 

z8d  =  z8  +  d8z 
c 

c   find  rotational  angles  of  element, 

c  based  on  the  vector  from  centroid  to  each  node  projected  into 
c   three  planes:  xy,  xz,  and  zy 
c   the  coordinates  of  the  centroid  of  element  =  (xcd,ycd, zed) 

xcdl  =  (Xl+x2+x3+x4+x5+x6+x7+x8) /8.d0 

ycdl  =  (yl+y2+y3+y4+y5+y6+y7+y8) /8.d0 

zcdl  =  (Zl+z2+z3+z4+z5+z6+z7+z8) /8.d0 
c 

xcd2  =  (Xld+x2d+x3d+x4d+x5d+x6d+x7d+x8d) /8.d0 
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ycd2  =  (yld+y2d+y3d+y4d+y5d+y6d+y7d+y8d) /8.d0 

zcd2  =  (Zld+z2d+z3d+z4d+z5d+z6d+z7d+z8d) /8.d0 
c 

call  thetafind  (xcdl,ycdl, zcdl,xcd2,ycd2, zcd2, 
+  xl , yl , zl , xld, yld, zld, ttdxl , ttdyl , ttdzl ) 

call  thetafind  (xcdl,ycdl, zcdl , xcd2 , ycd2 ,  zcd2 , 
+  x2 , y2 , z2 , x2d, y2d, z2d, ttdx2 , ttdy2 , ttdz2 ) 

call  thetafind  (xcdl,ycdl, zcdl,xcd2 ,ycd2,  zcd2, 
+  x3 ,y3 , z3 , x3d,y3d, z3d, ttdx3 ,  ttdy3 ,  ttdz3 ) 

call  thetafind  (xcdl,ycdl, zcdl,xcd2 ,ycd2,  zcd2, 
+  x4,y4, z4,x4d,y4d, z4d, ttdx4,  ttdy4,  ttdz4) 

call  thetafind  (xcdl,ycdl, zcdl,xcd2 ,ycd2,  zcd2, 
+  x5,y5, z5,x5d,y5d, z5d, ttdx5,  ttdy5,  ttdz5) 

call  thetafind  (xcdl,ycdl, zcdl,xcd2 ,ycd2, zcd2, 
+  x6,y6, z6,x6d,y6d, z6d, ttdx6,  ttdy6, ttdz6) 

call  thetafind  (xcdl,ycdl, zcdl,xcd2,ycd2,  zcd2, 
+  x7 , y7 , z7 ,x7d,y7d, z7d, ttdx7 ,  ttdy7  ,  ttdz7 ) 

call  thetafind  (xcdl,ycdl, zcdl , xcd2 , ycd2 ,  zcd2, 
+  x8 , y8 , z8 , x8d, y8d, z8d, ttdx8 , ttdy8 , ttdz8 ) 

c 

tx  =  (ttdxl+ttdx2+ttdx3+ttdx4+ttdx5+ttdx6+ttdx7+ttdx8) /8.d0 

ty  =  (ttdyl+ttdy2+ttdy3+ttdy4+ttdy5+ttdy6+ttdy7+ttdy8) /8.d0 

tz  =  (ttdzl+ttdz2+ttdz3+ttdz4+ttdz5+ttdz6+ttdz7+ttdz8) /8.d0 
c 

if  (dabs(tx)  .le.  1.0d-05)  tx  =  0 . OdO 

if  (dabs(ty)  .le.  1.0d-05)  ty  =  0 . OdO 

if  (dabs(tz)  .le.  1.0d-05)  tz  =  O.OdO 
c 

c   create  transformation  matrix  global  to  local 

c 

rll  =   dcos (tz) *dcos (ty)  -  dsin(tx) *dsin(ty) *dsin(tz) 

rml  =   dcos (ty) *dsin(tz)  +  dcos (tz) *dsin(ty) *dsin(tx) 

rnl  =  -dcos (tx) *dsin(ty) 

rl2  =  -dcos (tx) *dsin(tz) 

rm2  =   dcos (tz) *dcos (tx) 

rn2  =   dsin(tx) 

rl3  =   dsin(ty) *dcos (tz)  +  dcos (ty) *dsin(tx) *dsin(tz) 

rm3  =   dsin(ty) *dsin(tz)  -  dcos (ty) *dcos (tz) *dsin(tx) 

rn3  =   dcos (ty) *dcos (tx) 
c 

c   compute  local  node  coord,  and  displ.  (  rotation  from  global  to  local 
) 
c 

xll  =  rll*xl  +  rml*yl  +  rnl*zl 

yll  =  rl2*xl  +  rm2*yl  +  rn2*zl 

zll  =  rl3*xl  +  rm3*yl  +  rn3*zl 

xl2  =  rll*x2  +  rml*y2  +  rnl*z2 

yl2  =  rl2*x2  +  rm2*y2  +  rn2*z2 

zl2  =  rl3*x2  +  rm3*y2  +  rn3*z2 

xl3  =  rll*x3  +  rml*y3  +  rnl*z3 

yl3  =  rl2*x3  +  rm2*y3  +  rn2*z3 

zl3  =  rl3*x3  +  rm3*y3  +  rn3*z3 

xl4  =  rll*x4  +  rml*y4  +  rnl*z4 

yl4  =  rl2*x4  +  rm2*y4  +  rn2*z4 

zl4  =  rl3*x4  +  rm3*y4  +  rn3*z4 

xl5  =  rll*x5  +  rml*y5  +  rnl*z5 

yl5  =  rl2*x5  +  rm2*y5  +  rn2*z5 

zl5  =  rl3*x5  +  rm3*y5  +  rn3*z5 

xl6  =  rll*x6  +  rml*y6  +  rnl*z6 

yl6  =  rl2*x6  +  rm2*y6  +  rn2*z6 
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zl6  =  rl3 

*x6  +  rm3*y6 

+  rn3*z6 

xl7  =  rll 

*x7  +  rml*y7 

+  rnl*z7 

yl7  =  rl2 

*x7  +  rm2*y7 

+  rn2*z7 

zl7  =  rl3 

*x7  +  rm3  *y7 

+  rn3*z7 

xl8  =  rll 

*x8  +  rml*y8 

+  rnl*z8 

yl8  =  rl2 

*x8  +  rm2*y8 

+  rn2*z8 

zl8  =  rl3 

*x8  +  rm3*y8 

+  rn3*z8 

c 

disl(l) 

=  rll*dlx 

+ 

rml*dly  + 

rnl*dlz 

disl(2) 

=  rl2*dlx 

+ 

rm2*dly  + 

rn2*dlz 

disl(3) 

=   rl3*dlx 

+ 

rm3*dly  + 

rn3*dlz 

c 

disl(4) 

=  rll*d2x 

+ 

rml*d2y  + 

rnl*d2z 

disl(5) 

=   rl2*d2x 

+ 

rm2*d2y  + 

rn2*d2z 

disl(6) 

=   rl3*d2x 

+ 

rm3*d2y  + 

rn3*d2z 

c 

disl(7) 

=   rll*d3x 

+ 

rml*d3y  + 

rnl*d3z 

disl(8) 

=   rl2*d3x 

+ 

rm2*d3y  + 

rn2*d3z 

disl(9) 

=   rl3*d3x 

+ 

rm3*d3y  + 

rn3*d3z 

c 

disl(10) 

=   rll*d4x 

+ 

rml*d4y  + 

rnl*d4z 

disl(ll) 

=   rl2*d4x 

+ 

rm2*d4y  + 

rn2*d4z 

disl(12) 

=   rl3*d4x 

+ 

rm3*d4y  + 

rn3*d4z 

c 

disl(13) 

=  rll*d5x 

+ 

rml*d5y  + 

rnl*d5z 

disl(14) 

=   rl2*d5x 

+ 

rm2*d5y  + 

rn2*d5z 

disl(15) 

=  rl3*d5x 

+ 

rm3*d5y  + 

rn3*d5z 

c 

disl(16) 

=   rll*d6x 

+ 

rml*d6y  + 

rnl*d6z 

disl(17) 

=   rl2*d6x 

+ 

rm2*d6y  + 

rn2*d6z 

disl(18) 

=   rl3*d6x 

+ 

rm3*d6y  + 

rn3*d6z 

c 

disl(19) 

=  rll*d7x 

+ 

rml*d7y  + 

rnl*d7z 

disl(20) 

=   rl2*d7x 

+ 

rm2*d7y  + 

rn2*d7z 

disl(21) 

=   rl3*d7x 

+ 

rm3*d7y  + 

rn3*d7z 

c 

disl(22) 

=   rll*d8x 

+ 

rml*d8y  + 

rnl*d8z 

disl(23) 

=   rl2*d8x 

+ 

rm2*d8y  + 

rn2*d8z 

disl(24) 

=   rl3*d8x 

+ 

rm3*d8y  + 

rn3*d8z 

c 

c 

Initialized 

the  stresses 

and  internal  forces 

c 

do  100  j= 
epsx(j) 
epsy ( j ) 
epsz(j) 
epsxy(j) 
epsyz ( j ) 
epszx( j ) 
100  continue 

1,8 
=  O.dO 
=  O.dO 
=  O.dO 
=  O.dO 
=  O.dO 
=  O.dO 

c 

do  109  ijj  -  1,24 

f(ijj)  = 

O.dO 

109  continue 

at  each  Gauss ' s  point 


c 

c   coeff.  of  8  point  Gauss  Quadrature:  s (i) , t (i) ,r (i) ,w(i) ; i=l, 8 

c 

const  =  I.d0/dsqrt(3.0d0) 

s(l)  =  -const 

t(l)  =  -const 
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r(l)  =  const 


s(2 

= 

const 

t(2 

= 

-const 

r(2 

= 

const 

s(3 

= 

const 

t(3 

= 

const 

r(3 

= 

const 

s(4 

= 

-const 

t(4 

= 

const 

r(4' 

= 

const 

s(5] 

= 

-const 

t(5) 

= 

-const 

r(5] 

= 

-const 

s(6) 

= 

const 

t(6) 

= 

-const 

r(6) 

= 

-const 

s{7) 

= 

const 

t(7) 

= 

const 

r(7) 

= 

-const 

s(8] 

= 

-const 

t(8) 

= 

const 

r(8) 

= 

-const 

c 
c 
c 


do  199  iii  =  1,8 
w(iii)  =  l.dO 
.99  continue 

set  shape  function  Ni  =>  k  =  gauss ' spoint  no. 


do  201 
shpfl 
shpf2 
shpf3 
shpf4 
shpf5 
shpf  6 
shpf7 
shpf  8 


k= 


1,8 
(l.dO 


s(k)) 
d0+s(k) ) 
d0+s(k) ) 
dO-s(k) ) 
dO-s(k) ) 
d0+s(k) ) 
(l.d0+s(k) ) 
(l.dO-s(k) ) 


(1 

(1 
(1 
(1 
(1 


(1 
(1 
(1 
(1 
(1 
(1 
(1 
(1 


.dO-t(k) )*(1 
,d0-t(k) )*(1 
,dO+t(k) )* (1 
,dO+t(k) )*(1 
.dO-t(k) )* (1 
.dO-t(k) )* (1 
,dO+t(k) )* (1 
,dO+t(k) )* (1 


.d0+r(k) 
.d0+r(k) 
.d0+r(k) 
.d0+r(k) 
.dO-r(k) 
.dO-r(k) 
.dO-r(k) 
.dO-r(k) 


)/8.d0 
)/8.d0 
)/8.d0 
)/8.d0 
)/8.d0 
)/8.d0 
)/8.d0 
)/8.d0 


c 
c 
c 


derivatives  of  shape  function  w.r.t.  s,t,r 


sis 
s2s 
s3s 
s4s 
s5s 
s6s 
s7s 
s8s 


(l.dO-t(k)) 
(l.dO-t(k)) 
(l.d0+t(k) ) 
(l.d0+t(k) ) 
(l.dO-t(k)) 
(l.dO-t(k)) 
(l.dO+t(k)) 
(l.d0+t(k) ) 


(l.d0+r(k) ) 
(l.d0+r(k) ) 
(l.d0+r(k) ) 
(l.d0+r(k) ) 
(l.dO-r(k)) 
(l.dO-r(k)) 
(l.dO-r(k)) 
(l.dO-r(k)) 


.d0 
,d0 
,d0 
,d0 
.d0 
,d0 
,d0 
.d0 


sit  =  -(l.dO-s(k) )* (l.d0+r(k) )  /  8.d0 
s2t  =  -(l.d0+s(k) )*(l.d0+r(k) )  /  8.d0 
s3t  =   (l.d0+s(k) ) *(l.d0+r(k) )  /  8.d0 
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s4t  = 

(l.dO-s(k) )*(l.dO+r(k) ) 

/ 

8.  dO 

s5t  = 

-(l.dO-s(k) )*(l.dO-r(k)) 

/ 

8.  dO 

s6t  = 

-(l.dO+s(k) )*(l.dO-r(k) ) 

/ 

8.d0 

s7t  = 

(l.dO+s(k) )*(l.dO-r(k) ) 

/ 

8.d0 

s8t  = 

(l.dO-s(k) )*(l.dO-r(k) ) 

/ 

8.  dO 

c 

sir  = 

(l.dO-s(k) )*(l.dO-t(k) ) 

/ 

8.  dO 

s2r  = 

(l.dO+s(k) )*(l.dO-t(k)) 

/ 

8.d0 

s3r  = 

(l.dO+s(k) )*(l.dO+t(k) ) 

/ 

8.  dO 

s4r  = 

(l.dO-s(k) )*(l.dO+t(k)) 

/ 

8.d0 

s5r  = 

-(l.dO-s(k) )*(l.dO-t(k)) 

/ 

8.d0 

s6r  = 

-(l.dO+s(k) )*(l.dO-t(k)) 

/ 

8.d0 

s7r  = 

-(l.dO+s(k) )*(l.dO+t(k)) 

/ 

8.  dO 

s8r  = 

-(l.dO-s(k) )*(l.dO+t(k) ) 

/ 

8.  dO 

c 

c  Jacobian  Matrix  Element ( Jik) 

c 

xJll  = 

Sls*xll+s2s*xl2+s3s*xl3 

^s4s*xl4+ 

+ 

S5s*xl5+s6s*xl6+s7s*xl7 

*s8s*xl8 

xJ12  = 

sls*yll+s2s*yl2+s3s*yl3 

i-s4s*yl4+ 

+ 

S5s*yl5+s6s*yl6+s7s*yl7 

*-s8s*yl8 

xJ13  = 

Sls*zll+s2s*zl2+s3s*zl3 

4-s4s*zl4+ 

+ 

S5s*zl5+s6s*zl6+s7s*zl7 

*s8s*zl8 

xJ21  = 

Slt*xll+s2t*xl2+s3t*xl3 

fs4t*xl4+ 

+ 

S5t*xl5+s6t*xl6+s7t*xl7 

t-s8t*xl8 

xJ22  = 

slt*yll+s2t*yl2+s3t*yl3 

<-s4t*yl4  + 

+ 

S5t*yl5+s6t*yl6+s7t*yl7+s8t*yl8 

xJ23  = 

Slt*zll+s2t*zl2+s3t*zl3 

fs4t*zl4+ 

+ 

S5t*zl5+s6t*zl6+s7t*zl7 

4-s8t*zl8 

xJ31  = 

Slr*xll+s2r*xl2+s3r*xl3 

fs4 

ir*xl4+ 

+ 

s5r*xl5+s6r*xl6+s7r*xl7+s£ 

Ir*xl8 

xJ32  = 

Slr*yll+s2r*yl2+s3r*yl3 

fs4 

Lr*yl4+ 

+ 

S5r*yl5+s6r*yl6+s7r*yl7 

fsE 

!r*yl8 

xJ33  = 

Slr*zll+s2r*zl2+s3r*zl3 

t-s4 

ir*zl4+ 

+ 

S5r*zl5+s6r*zl6+s7r*zl7 

(-sE 

r*zl8 

c 
c 


Jacobian  =  det | J | 


xJ12*xJ23*xJ31 
xJ32*xJ23*xJll 


Node  No . 
Node  No . 


down     =  (xJll*xJ22*xJ33  + 

up       =  (xJ31*xJ22*xJl3  + 

det  J     =  down  -  up 
check  negative  jacobian 
if  (detJ  .le.  0.0)  then 
if  (k  .eq.  1)  then 

print  *,   ' Negative  Jacobian, 

write (6,*)' Negative  Jacobian, 

else  if  (k  .eq.  2)  then 

print  *,   ' Negative  Jacobian,  Node  No. 

write (6,*)' Negative  Jacobian,  Node  No. 

else  if  (k  .eq.  3)  then 

print  *,   ' Negative  Jacobian,  Node  No. 

write (6,*)' Negative  Jacobian,  Node  No. 

else  if  (k  .eq.  4)  then 

print  *,   ' Negative  Jacobian,  Node  No 

write (6,*)' Negative  Jacobian,  Node  No. 

else  if  (k  .eq.  5)  then 

print  *,   ' Negative  Jacobian,  Node  No. 

write(6,*)' Negative  Jacobian,  Node  No. 

else  if  (k  .eq.  6)  then 


+  xJ13*xJ21*xJ32) 
+  xJ33*xJ21*xJ12) 


nl  , 
nl  , 

n2  , 
n2  , 

n3  , 
n3  , 

,  n4 
n4  , 

n5  , 
n5  , 


J= 
J= 

J= 
J= 

J= 
J= 

1  J= 


J= 
J= 


detJ 
det  J 

det  J 
detJ 

det  J 
det  J 

,  detJ 
detJ 

det  J 
det  J 
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c 
c 
c 


c 
c 

c 


print  * 
write (6 

else  if 
print  * 
write (6 

else  if 
print  * 
write (6 

end  if 

mode  (12,  i) 
go  to  13 
endif 


' Negative  Jacobian,  Node  No.  = 

)' Negative  Jacobian,  Node  No.  = 

(k  .eq.  7)  then 

' Negative  Jacobian,  Node  No.  = 

) ' Negative  Jacobian,  Node  No.  = 

(k  .eq.  8)  then 

■ Negative  Jacobian,  Node  No.  = 

Negative  Jacobian,  Node  No.  = 


)  '-- 


'  ,  n6  , 

J=  ' 

,  detJ 

1  ,  n6  , 

J=  ' 

,  detJ 

\  n7  , 

J=  ' 

,  detJ 

'  ,  n7  , 

J=  ' 

,  detJ 

■  ,  n8  , 

J=  ' 

,  detJ 

'  ,  n8  , 

J=  ' 

,  detJ 

11. dO 


form  inverse  jacobian 


all 
al2 
al3 
a21 
a22 
a23 
a31 
a32 
a33 


(xJ22*xJ33 
(xJ12*xJ33 
(xJ12*xJ23 
(xJ21*xJ33 
(xJll*xJ33 
(xJll*xJ23 
(xJ21*xJ32 
(xJll*xJ23 
(xJll*xJ22 


xJ32*xJ23) /detJ 
xJ32*xJ13) /detJ 
xJ22*xJ13) /detJ 
xJ31*xJ23) /detJ 
xJ31*xJ13) /detJ 
xJ21*xJ13) /detJ 
xJ31*xJ22) /detJ 
xJ21*xJ13) /detJ 
xJ21*xJ12) /detJ 


form  B  matrix 


blx  =  all*sls+al2*slt+al3*slr 
b2x  =  all*s2s+al2*s2t+al3*s2r 
b3x  =  all*s3s+al2*s3t+al3*s3r 
b4x  =  all*s4s+al2*s4t+al3*s4r 
b5x  =  all*s5s+al2*s5t+al3*s5r 
b6x  =  all*s6s+al2*s6t+al3*s6r 
b7x  =  all*s7s+al2*s7t+al3*s7r 
b8x  =  all*s8s+al2*s8t+al3*s8r 

bly  =  a21*sls+a22*slt+a23*slr 
b2y  =  a21*s2s+a22*s2t+a23*s2r 
b3y  =  a21*s3s+a22*s3t+a23*s3r 
b4y  =  a21*s4s+a22*s4t+a23*s4r 
b5y  =  a21*s5s+a22*s5t+a23*s5r 
b6y  =  a21*s6s+a22*s6t+a23*s6r 
b7y  =  a21*s7s+a22*s7t+a23*s7r 
b8y  =  a21*s8s+a22*s8t+a23*s8r 

biz  =  a31*sls+a32*slt+a33*slr 
b2z  =  a31*s2s+a32*s2t+a33*s2r 
b3z  =  a31*s3s+a32*s3t+a33*s3r 
b4z  =  a31*s4s+a32*s4t+a33*s4r 
b5z  =  a31*s5s+a32*s5t+a33*s5r 
b6z  =  a31*s6s+a32*s6t+a33*s6r 
b7z  =  a31*s7s+a32*s7t+a33*s7r 
b8z  =  a31*s8s+a32*s8t+a33*s8r 


do  21  ie=l, 6 
do  22  le=l,24 
b(ie,le)  =  O.dO 

22  continue 

21  continue 

b(l,l)   =  blx 
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b(l 
b(l 
b(l 
b(l 
b(l 
b(l 
b(l 
b(2 
b(2 
b(2 
b(2 
b(2 
b(2 
b(2 
b(2 
b(3 
b(3 
b(3 
b(3 
b(3 
b(3 
b(3 
b(3 
b(4 
b(4 
b(4 
b(4 
b(4 
b(4 
b(4 
b(4 
b(4 
b(4 
b(4 
b(4 
b(4 
b(4 
b(4 
b(4 
b(5 
b(5 
b(5 
b(5 
b(5 
b(5 
b(5 
b(5 
b(5 
b(5 
b(5 
b(5 
b(5 
b(5 
b(5 
b(5 
b(6 
b(6 
b(6 
b(6 
b(6 


4) 

=  b2x 

7) 

=  b3x 

10) 

=  b4x 

13) 

=  b5x 

16) 

=  b6x 

19) 

=  b7x 

22) 

=  b8x 

2) 

=  bly 

5) 

=  b2y 

8) 

=  b3y 

11) 

=  b4y 

14) 

=  b5y 

17) 

=  b6y 

20) 

=  b7y 

23) 

=  b8y 

3) 

=  biz 

6) 

=  b2z 

9) 

=  b3z 

12) 

=  b4z 

15) 

=  b5z 

18) 

=  b6z 

21) 

=  b7z 

24) 

=  b8z 

1) 

=  bly 

2) 

=  blx 

4) 

=  b2y 

5) 

=  b2x 

7) 

=  b3y 

8) 

=  b3x 

10) 

=  b4y 

11) 

=  b4x 

13) 

=  b5y 

14) 

=  b5x 

16) 

=  b6y 

17) 

=  b6x 

19) 

=  b7y 

20) 

=  b7x 

22) 

=  b8y 

23) 

=  b8x 

2) 

=  biz 

3) 

=  bly 

5) 

=  b2z 

6) 

=  b2y 

8) 

=  b3z 

9) 

=  b3y 

11) 

=  b4z 

12) 

=  b4y 

14) 

=  b5z 

15) 

=  b5y 

17) 

=  b6z 

18) 

=  b6y 

20) 

=  b7z 

21) 

=  b7y 

23) 

=  b8z 

24) 

=  b8y 

1) 

=  biz 

3) 

=  blx 

4) 

=  b2z 

6) 

=  b2x 

7) 

=  b3z 
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b(6,9)  =  b3x 

b(6,10)  =  b4z 

b(6,12)  =  b4x 

b(6,13)  =  b5z 

b(6,15)  =  b5x 

b(6,16)  =  b6z 

b(6,18)  =  b6x 

b(6,19)  =  b7z 

b(6,21)  =  b7x 

b(6,22)  =  b8z 

b(6,24)  =  b8x 


c   compute  local  element  strains  of  each  Gauss ' s  point 
c 

do  300  j=l,24 
epsx(k)  =  epsx(k)  +  b(l, j ) *disl ( j ) 
epsy(k)  =  epsy(k)  +  b(2, j ) *disl ( j ) 
epsz(k)  =  epsz(k)  +  b(3, j ) *disl ( j ) 
epsxy(k)  =  epsxy(k)  +  b(4, j ) *disl ( j ) 
epsyz(k)  =  epsyz(k)  +  b(5, j ) *disl ( j ) 
epszx(k)  =  epszx(k)  +  b(6, j ) *disl ( j ) 
300  continue 
c 

c  assign  value  of  new  epslon  (epslonN(..) 
c 

m  =  ( (i-l)*8+k-l)*6 
epslonN(m+l)  =  epsx(k) 
epslonN(m+2)  =  epsy(k) 
epslonN (m+3)  =  epsz (k) 
epslonN(m+4)  =  epsxy(k) 
epslonN (m+5)  =  epsyz(k) 
epslonN (m+6)  =  epszx(k) 
c 

c  calculate  delta  epslon 
c 

do  352  in  =  1,6 
Depslon(in)  =  epslonN (m+in)  -  epslonP(m+in) 
352  continue 
c 

c  material  properties  for  element 
c 

mtyp    =  int  (mode  (9,  i)  ) 

mtyp2    =  int (e(l, mtyp) ) 

eyng    =  e (3, mtyp) 

if  (eyng  .le.  O.OdO)  goto  13 

poisson  =  e (4, mtyp) 

ft      =  e(5,mtyp) 

coh     =  e(6,mtyp) 

phiangle  =  e (7,mtyp) *phi/180 .d0 

kickPL   =  int(e(8,mtyp) ) 

eyngt    =  e(9,mtyp) 

beta     =  e( 10, mtyp) 

tau      =  e(ll,l) 


c 
c 

c 

c 
c 
c 


call  constitutive  coefficient  from  subroutine, 
call  elastd3d  (eyng, poisson, ed) 
calculate  trial  stress (sigma  =  stresses  at  gauss  pts) 
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353 


ny  =  (i-1) *8  +  k 

PLr  =  PLrP(ny) 

do  353  m=l,6 
ny  =  ( (i-l)*8+k-l)*6+m 
sigma(m)  =  sigmaP(ny) 
PLalpha(m)  =  PLalphaP(ny) 

continue 


if  (kickPL  .eq.  1)  then 
do  402  m  =  1,6 

do  502  n  =  1,6 

ed(m,n)  =  exp(-time/tau) *ed(m,n) 
502  continue 
402  continue 

end  if 

do  400  m  =  1, 6 

do  500  n  =  1, 6 

sigma(m)  =  sigma (m)  +  ed(m,n) *Depslon(n) 
500    continue 
400  continue 


c 
c 

c 
c 
c 


c 
c 
c 
c 
c 


sigmean  =  sigma (1) +sigma(2) +sigma (3) /3 . OdO 


apply  plasticity  model  :  Drucker-Prager  yield  criterion 

if(kickPL  .gt.  1) 
+call  PLmodel (sigma, eyng,poisson, ft, coh,phiangle, kickPL, 
+  eyngt,beta, PLalpha, PLr, sigmean , mtyp2 , elplas, i,k) 

transform  local  nodal  stresses  into  global 

— sig(i,k)  =  sigma{i=x- (1) ,y- (2) , z- (3) ,xy- (4) , 

yz- (5) ,xz- (6) } , gauss  pts(k=l,8) 


te(l,l 
te(l,2 
te(l,3 
ted, 4 
ted,  5 
ted, 6 
te(2,l 
te(2,2 
ted,  3 
ted, 4 
ted, 5 
ted,  6 
te(3,l 
ted,  2 
te(3, 3 
ted,  4 
te(3, 5 
ted, 6 
te(4,l 
te(4,2 
te(4,3 
te(4,4 
te(4,5 
te(4, 6 


=rll*rll 

=rml*rml 

=rnl*rnl 

=rll*rml 

=rml*rnl 

=rnl*rll 

=rl2*rl2 

=rm2*rm2 

=rn2*m2 

=rl2*rm2 

=rm2*rn2 

=rn2*rl2 

=rl3*rl3 

=rm3*rm3 

=rn3*rn3 

=rl3*rm3 

=rm3*rn3 

=rn3*rl3 

=2.0d0*rll*rl2 

=2.0d0*rml*rm2 

=2.0d0*rnl*rn2 

=rll*rm2+rl2*rml 

=rml * rn2 +rm2 * rnl 

=rnl*rl2+rn2*rll 
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te(5,l)=2.0d0*rl2*rl3 
te ( 5 ,  2 ) =2 . OdO *rm2  *rm3 
te(5,3)=2.0dO*m2*rn3 
te (5, 4) =rl2*rm3+rl3*rm2 
te(5, 5)=rm2*rn3+rm3*rn2 
te ( 5 , 6 ) =rn2*rl3+rn3*rl2 
te(6,l)=2.0d0*rl3*rll 
te  ( 6 ,  2 )  =2  .  OdO *rm3  *rml 
te ( 6 , 3 ) =2 . OdO  *rn3  *rnl 
te(6,4)=rl3*rml+rll*rm3 
te (6, 5) =rm3*rnl+rml*rn3 
te (6, 6) =rn3*rll+rnl*rl3 
c 

do  777  mi=l, 6 
sig(k,mi)  =  O.OdO 
do  810  ni=l,6 

sig(k,mi)  =  sig(k,mi)  +  te (ni,mi) *sigma(ni) 
810  continue 
777  continue 
c 

c  compute  local  element  internal  nodal  forces 
c 

do  600  111=1,24 
do  700  n=l,6 
f (m)  =  f(m)  +  w(k) *b(n,m) *sigma(n) *detJ 
700  continue 
600  continue 
c 

PLrN( (i-l)*8+k)  =  PLr 
do  800  m=l, 6 
ny  =  ( (i-l)*8+k-l)*6+m 
sigmaN(ny)  =  sigma(m) 
PLalphaN(ny)  =  PLalpha(m) 
800  continue 
201  continue 
c 

c  Extrapolate  stress  from  Gauss's  pts  to  node 
c  and  calculate  average  nodal  stress 
c 

call  stress (sig,nl,n2,n3,n4,n5,n6,n7,n8, si, s2 , s3, s4, s5, s6, s7) 

c 

c   transform  local  internal  nodal  forces  to  global  internal  nodal  forces 

c 

detT  =  rll*rm2*rn3+rl2*rm3*rnl+rl3*rml*rn2 
+       -rl3*rm2*rnl-rl2*rml*rn3-rll*rm3*rn2 

sll  =   (rm2*rn3  -  rm3*rn2) /detT 

sml  =  -(rml*rn3  -  rm3*rnl) /detT 

snl  =   (rml*rn2  -  rm2*rnl) /detT 

sl2  =  -(rl2*rn3  -  rl3*rn2) /detT 

sm2  =   (rll*rn3  -  rl3*rnl) /detT 

sn2  =  -(rll*rn2  -  rl2*rnl) /detT 

sl3  =   (rl2*rm3  -  rl3*rm2) /detT 

sm3  =  -(rll*rm3  -  rl3*rml) /detT 

sn3  =   (rll*rm2  -  rl2*rml) /detT 
c 

fix  =  sll*f (1)   +  sml*f(2)   +  snl*f (3) 

fly  =  sl2*f (1)   +  sm2*f(2)   +  sn2*f (3) 

flz  =  sl3*f(l)   +  sm3*f(2)   +  sn3*f(3) 

f2x  =  sll*f (4)   +  sml*f (5)   +  snl*f(6) 
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f2y 

= 

sl2*f (4) 

+ 

sm2*f (5) 

+ 

sn2*f (6) 

f2z 

= 

sl3*f  (4) 

+ 

sm3*f (5) 

+ 

sn3*f (6) 

f3x 

= 

sll*f (7) 

+ 

sml*f (8) 

+ 

snl*f (9) 

f3y 

= 

sl2*f (7) 

+ 

sm2*f (8) 

+ 

sn2*f (9) 

f3z 

= 

sl3*f (7) 

+ 

sm3*f (8) 

+ 

sn3*f (9) 

f4x 

= 

sll*f (10) 

+ 

sml*f (11 

+ 

snl*f (12) 

f4y 

= 

sl2*f (10) 

+ 

sm2*f (11 

+ 

sn2*f (12) 

f4z 

= 

sl3*f (10) 

+ 

sm3*f (11 

+ 

sn3*f (12) 

f5x 

= 

sll*f (13) 

+ 

sml*f (14 

+ 

snl*f (15) 

f5y 

= 

sl2*f (13) 

+ 

sm2*f (14 

+ 

sn2*f (15) 

f5z 

= 

sl3*f (13) 

+ 

sm3*f (14 

+ 

sn3*f (15) 

f6x 

= 

sll*f (16) 

+ 

sml*f (17 

+ 

snl*f (18) 

f6y 

= 

sl2*f (16) 

+ 

sm2*f (17 

+ 

sn2*f (18) 

f6z 

= 

sl3*f (16) 

+ 

sm3*f (17 

+ 

sn3*f (18) 

fix 

= 

sll*f (19) 

+ 

sml*f (20 

+ 

snl*f (21) 

fly 

= 

sl2*f (19) 

+ 

sm2*f (20 

+ 

sn2*f (21) 

f7z 

= 

sl3*f (19) 

+ 

sm3*f (20 

+ 

sn3*f (21) 

f8x 

= 

sll*f (22) 

+ 

sml*f (23 

+ 

snl*f (24) 

f8y 

= 

sl2*f (22) 

+ 

sm2*f (23 

+ 

sn2*f (24) 

f8z 

= 

sl3*f (22) 

+ 

sm3*f (23 

+ 

sn3*f (24) 

c 
c 

assemble  to  global  nodal  : 

force  m 

c 

pint (ml  +  1 

= 

pint (ml+1 

+ 

fix 

pint (ml+2 

= 

pint (ml+2 

+ 

fly 

pint (ml+3 

= 

pint (ml+3 

+ 

flz 

pint (m2+l 

= 

pint (m2+l 

+ 

f2x 

pint (m2+2 

= 

pint (m2+2 

+ 

f2y 

pint (m2+3 

= 

pint (m2+3 

+ 

f2z 

pint (m3  +  l 

= 

pint (m3+l 

+ 

f3x 

pint (m3+2) 

= 

pint (m3+2 

+ 

f3y 

pint (m3+3) 

= 

pint (m3+3 

+ 

f3z 

pint (m4  +  l 

= 

pint (m4+l 

+ 

f4x 

pint (m4+2 

= 

pint (m4+2 

+ 

f4y 

pint (m4+3 

= 

pint (m4+3 

+ 

f4z 

pint (m5+l) 

= 

pint (m5+l 

+ 

f5x 

pint (m5+2 

= 

pint (m5+2 

+ 

f5y 

pint (m5+3 

= 

pint (m5+3 

+ 

f5z 

pint (m6  +  l 

= 

pint (m6+l 

+ 

f6x 

pint (m6+2 

= 

pint (m6+2 

+ 

f6y 

pint (m6+3 

= 

pint (m6+3 

+ 

f6z 

pint (m7  +  l 

= 

pint (m7+l 

+ 

f7x 

pint (m7+2 

= 

pint (m7+2 

+ 

f7y 

pint (m7  +  3 

= 

pint (m7+3 

+ 

f7z 

pint (mS+l1 

= 

pint (m8+l 

+ 

f8x 

pint (m8+21 

= 

pint (m8+2 

+ 

f8y 

pint (m8+3' 

= 

pint (m8+3 

+ 

f8z 

13  return 

end 

pint  (internal) 


it*********************************************************************** 

*  Drucker-Prager  yield  criterion  with  non-associate  flow  rule  and      * 

*  linear  combination  of  isotropic  and  kinematic  hardening/softening    * 

*  Krieq  and  Key's  radial-return  algorithm  for  elastoplastic  case       * 

*  (  kickPL  =  3  for  Drucker-Prager)  * 

*  (  beta  =  0.  for  kinematic  &  =1  for  isotropic  hardening)      * 
a*********************************************************************** 


subroutine  PLmodel (sigma, eyng,poisson, ft, coh,phiangle, kickPL, 
+  eyngt,beta, PLalpha, PLr, sigmean,mtyp2 , elplas, i,k) 
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implicit  real*8  (a-h,o-z) 

dimension  PLalpha ( 1 ) , sigma ( 1 ) , xi ( 6 ) , elplas ( 1 ) 
c 

Ep  =  eyngt/ (1. OdO-eyngt/eyng) 
gshear  =  eyng/ (2 . OdO* (1 . OdO+poisson) ) 
cl  =  2.0d0*gshear  +  Ep*2 . OdO/3 . OdO 
c2  =  2.0dO/3.0dO*beta*Ep 
c3  =  2.0dO/3.0dO*(1.0dO-beta)*Ep 
c4  =  2.0d0*gshear 
c 

c  calculate  xi-trial 
c 

do  10  j=l,6 
xi(j)  =  sigma(j)  -  PLalpha(j) 
10  continue 
c 

c  calculate  deviatoric  part  of  xi 
c 

ximean  =  (xi (1) +xi (2) +xi (3) ) /3 . OdO 
c 

c  calculate  the  radius  of  cross  section  of  cone/cylinder 
c 

if  ( (mtyp2  .eq.  2) .and. (kickPL  .eq.  2))  phiangle  =  0 . OdO 
if(kickPL  .eq.  2)  then 
b  =  O.OdO 

ak  =  ft/dsqrt (3. OdO) 
endif 

if  (kickPL  .eq.  3)  then 
b  =  dsqrt (12. dO) *dsin (phiangle) / (3 . OdO+dsin (phiangle) ) 
ak  =  dsqrt (12 .dO) * (-coh) *dcos (phiangle) / (3 .dO+dsin (phiangle) ) 
endif 

r  =  dabs (dsqrt (2.d0) * (ak+b*sigmean) ) 
if  (r  .gt.  PLr)  PLr  =  r 
c 

do  20  j  =1,3 
xi  ( j  )  =  xi  ( j  )  -  ximean 
20  continue 
c 

c  check  if  elastic 
c 

xxi  =  xi(l)*xi(l)+xi(2)*xi(2)+xi(3)*xi(3)  + 
+      2 . dO  * (xi ( 4 ) *xi ( 4 ) +xi ( 5 ) *xi ( 5 ) +xi ( 6 ) *xi ( 6 ) ) 
yn  =  PLr* PLr 

if  (xxi  .le.  yn)  go  to  100 
c 

c  Plastic  phase:  calculate  unit  normal--N  (also  store  in  "xi") 
c 

xxi  =  dsqrt (xxi) 
do  30  j  =1,6 
xi  ( j  )  =  xi  ( j  )  /xxi 
30  continue 
c 

c  calculate  lambda-tilde  ("A") 
c 

A  =  (xxi-PLr) /cl 
c 

c  update  stresses 
c 

PLr  =  PLr  +  c2*A 
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tl  =  c3*A 
t2  =  c4*A 
do  40  j=l,6 
PLalpha(j)  =  PLalpha(j)  +  tl*xi(j) 
sigma(j)  =  sigma(j)  -  t2*xi(j) 
40  continue 
c 
c  record  the  element  number  if  plastic 


c 


j  =  8*(i-l)+k 

elplas(j)  =  dble(i)  +  0 . Id0*dble(k) 


100  return 
end 

c 

************************************************************************ 

*  * 

*  compute  element  internal  nodal  streses (sigma)  by   * 

*  8-node  solid  isoparametric  element  (8-node  hexahedral  element)       * 

*  * 

************************************************************************ 

c 

subroutine  stress (sig,nl,n2,n3,n4,n5,n6,n7,n8, 
+  si, s2 , s3, s4, s5, s6, s7) 

c 

implicit  real*8  (a-h,o-z) 

dimension  x(8) ,y(8) ,z(8) 

dimension  sig ( 8 , 6 ) , signodeX ( 8 ) , signodeY ( 8 ) , signodeZ ( 8 ) 

dimension  signodeXY ( 8 ) , signodeYZ ( 8 )  , signodeXZ ( 8 ) 

dimension  sN(8,  8)  ,  si  (1)  ,  s2  (1)  ,  s3  (1)  ,  s4  (1)  ,  s5  (1)  ,  s6  (1)  ,  s7  (1) 
c 

common  /box  1/  iprob, delta, alpha, toler, gravity 
c 

c  Interpolate  stresses  from  gauss  points  to  nodes 
c  (  si(idof,io)  ;  idof=dof,  io=  local  node  no.  ) 
c 

constant  =  dsqrt(3.0d0) 

x(l)  =  -constant 

y(l)  =  -constant 

z(l)  =   constant 
c 

x(2)  =  constant 

y(2)  =  -constant 

z  (2)  =   constant 
c 

x(3)  =   constant 

y(3)  =   constant 

z(3)  =   constant 
c 

x(4)  =  -constant 

y(4)  =   constant 

z(4)  =   constant 
c 

x(5)  =  -constant 

y(5)  =  -constant 

z(5)  =  -constant 
c 

x(6)  =   constant 

y(6)  =  -constant 
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z(6) 


-constant 


x(7) 
y(7) 

z(7) 

x(8) 
y(8) 
z(8) 


cons 

cons 

-cons 

-cons 

cons 

-cons 


tant 
tant 
tant 

tant 
tant 
tant 


203 


do  203  j 
sN(j,l)  = 
sN(j,2)  = 
sN(j,3)  = 
sN(j,4)  = 
sN(j,5)  = 
sN(j,6)  = 
sN(j,7)  = 
sN(j,8)  = 
continue 


=  1 

■    (1 

■■    (1 

>    (1 

(1 

(1 

(1 

(1 

(1 


,8 

.dO-x(j)) 

.d0+x(j)) 

-d0+x(j)) 

•dO-x(j)) 

.dO-x(j)) 

.d0+x(j)) 

■d0+x(j)) 

.dO-x(j)) 


(1 

(1 
(1 
(1 
(1 
(1 
(1 
(1 


dO-y(j) 
dO-y(j) 
d0+y(j) 
d0+y(j) 
dO-y(j) 
dO-y(j) 
d0+y(j) 
d0+y(j) 


(l.d0+z(j) )/8.d0 
(l.d0+z(j) )/8.d0 
(l.d0+z(j) )/8.d0 
(l.d0+z(j) )/8.d0 
(l.dO-z(j) )/8.d0 
(l.dO-z(j) )/8.d0 
(l.dO-z(j) )/8.d0 
(l.dO-z(j) )/8.d0 


502 
501 


do  501  j  =  1,8 
signodeX(j)  = 
signodeY(j)  = 
signodeZ(j)  = 
signodeXY(j)  = 
signodeYZ (j )  = 
signodeXZ(j)  = 
do  502  k  =  1, 
signodeX( j ) 
signodeY( j ) 
signodeZ ( j ) 
signodeXY( j 
signodeYZ ( j 
signodeXZ ( j 
continue 
continue 


sl(nl) 
sl(n2) 
sl(n3) 
sl(n4) 
sl(n5) 
sl(n6) 
sl(n7) 
sl(n8) 


sl(nl) 
sl(n2) 
sl(n3) 
sl(n4) 
sl(n5) 
sl(n6) 
sl(n7) 
sl(n8) 


O.OdO 

O.OdO 

O.OdO 
O.OdO 
O.OdO 
O.OdO 

8 
=  signodeX(j)  + 
=  signodeY(j)  + 
=  signodeZ (j)  + 

}  =  signodeXY(j) 

)  =  signodeYZ (j) 

)  =  signodeXZ (j) 


signodeX (1) 
signodeX(2) 
signodeX (3 ) 
signodeX (4) 
signodeX (5) 
signodeX (6) 
signodeX (7) 
signodeX(8) 


sN(j,k)*sig(k,l) 
sN(j,k) *sig(k,2) 
sN(j,k)*sig(k,3) 
+  sN(j,k)*sig(k,4) 
+  sN(j,k)*sig(k,5) 
+  sN(j,k)*sig(k,6) 


s2(nl) 
s2 (n2) 
s2(n3) 
s2(n4) 
s2(n5) 
s2(n6) 
s2(n7) 
s2(n8) 


s2(nl) 
s2(n2) 
s2(n3) 
s2(n4) 
s2(n5) 
s2(n6) 
s2(n7) 
s2(n8) 


signodeY(l) 
signodeY(2) 
signodeY(3) 
signodeY(4) 
signodeY(5) 
signodeY(6) 
signodeY(7) 
signodeY(8) 


s3(nl)  =  s3(nl)  +  signodeZ (1) 
s3(n2)  =  s3(n2)  +  signodeZ (2) 
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s3(n3) 

= 

s3(n3 

+ 

signodeZ (3) 

s3 (n4) 

= 

s3  (n4 

+ 

signodeZ (4) 

s3(n5) 

= 

s3  (n5 

+ 

signodeZ (5) 

s3 (n6) 

= 

s3  (n6 

+ 

signodeZ (6) 

s3(n7) 

= 

s3  (n7 

+ 

signodeZ (7) 

s3 (n8) 

= 

s3  (n8 

+ 

signodeZ (8) 

s4 (nl) 

s 

s4  (nl 

+ 

signodeXY(l) 

s4 (n2) 

= 

s4  (n2 

+ 

signodeXY(2) 

s4 (n3) 

= 

s4  (n3 

+ 

signodeXY(3) 

s4 (n4) 

= 

s4  (n4 

+ 

signodeXY(4) 

s4 (n5) 

= 

s4  (n5 

+ 

signodeXY(5) 

s4 (n6) 

= 

s4  (n6 

+ 

signodeXY(6) 

s4 (n7) 

= 

s4  (n7 

+ 

signodeXY(7) 

s4 (n8) 

= 

s4(n8 

+ 

signodeXY(8) 

s5(nl) 

= 

s5  (nl 

+ 

signodeYZ (1) 

s5 (n2) 

= 

s5  (n2 

+ 

signodeYZ (2) 

s5 (n3) 

= 

s5(n3 

+ 

signodeYZ (3) 

s5 (n4) 

= 

s5  (n4 

+ 

signodeYZ (4) 

s5(n5) 

= 

s5(n5 

+ 

signodeYZ (5) 

s5 (n6) 

= 

s5  (n6 

+ 

signodeYZ (6) 

s5 (n7) 

= 

s5(n7 

+ 

signodeYZ (7) 

s5(n8) 

= 

s5  (n8 

+ 

signodeYZ (8) 

s6 (nl) 

= 

s6  (nl 

+ 

signodeXZ (1) 

s6 (n2) 

= 

s6  (n2 

+ 

signodeXZ (2) 

s6 (n3) 

= 

s6  (n3 

+ 

signodeXZ (3) 

s6 (n4) 

= 

s6  (n4 

+ 

signodeXZ (4) 

s6 (n5) 

= 

s6  (n5 

+ 

signodeXZ (5) 

s6 (n6) 

= 

s6  (n6 

+ 

signodeXZ (6) 

s6 (n7) 

= 

s6  (n7 

+ 

signodeXZ (7) 

s6 (n8) 

= 

s6  (n8 

+ 

signodeXZ (8) 

s7 (nl) 

= 

s7(nl 

+ 

l.OdO 

s7 (n2) 

= 

s7  (n2 

+ 

l.OdO 

s7(n3) 

= 

s7(n3 

+ 

l.OdO 

s7 (n4) 

= 

s7  (n4 

+ 

l.OdO 

s7(n5) 

= 

s7(n5 

+ 

l.OdO 

s7 (n6) 

= 

s7  (n6 

+ 

l.OdO 

s7 (n7) 

= 

s7(n7 

+ 

l.OdO 

s7 (n8) 

= 

s7(n8 

+ 

l.OdO 

13  return 

end 
************************************************************************ 

*  * 

*  print  out  final  results  requested  at  specific  points  &  directions   * 

*  * 

************************************************************************ 

c 

subroutine  prout  (numout, rkout,d, v,a, si, s2, s3 , s4, s5, s6, s7,ndof , 
+  nstep, time,proutv) 


c   This  subroutine  controls  output  of  displacement,  velocity, 

c  acceleration,  and  stresses. 

c 

implicit  real*8  (a-h,o-z) 

dimension  d(ndof , 1)  , vfndof , 1) ,a (ndof , 1) 

dimension  rkout (3,1), proutv (I),sl(l),s2(l),s3(l),s4(l) 


dimension  s5 (1) , s6 (1) , s7 (1) 
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do  100  i2= 
node  =  int 
idva  =  int 
idof  =  int 
if  (idva 
if  (idva 
if  (idva 
if  (idva 
if  (idof 
if  (idof 
if  (idof 
if  (idof 
if  (idof 
if  (idof 
end  if 
100  continue 

write  (*, 
200  format ('  [ 
write(7,30 
300  format (lx, 
19x, 
19x, 
19x, 
19x, 
19x, 
19x, 
19x, 
19x, 
19x, 
19x, 
19x, 
19x, 
19x, 
19x, 


1 , numout 
(rkout(l,i2)) 
(rkout(2,i2)) 
(rkout(3,i2) ) 


.eq. 
.eq. 
.eq. 
.eq. 

.eq. 

.eq. 

.eq. 

.eq. 

.eq. 

.eq. 


0 )  proutv ( i 

1 )  proutv ( i 

2)  proutv (i 

3 )  then 

1 )  proutv ( 
proutv ( 
proutv ( 
proutv ( 
proutv ( 
proutv ( 


2)  =  d( idof, node) 
2)  =  v( idof, node) 
2)  =  a (idof, node) 


2) 
3) 
4) 
5) 
6) 


i2) 
i2) 
12) 
i2) 
i2) 
i2) 


si (node) 
s2 (node) 
s3 (node) 
s4 (node) 
s5 (node) 
s6 (node) 


/s7 (node) 
/s7 (node) 
/s7 (node) 
/s7 (node) 
/s7 (node) 
/s7 (node) 


200)  nstep 


+ 
+ 
+ 
+ 
+ 
+ 
+ 
+ 
+ 
+ 
+ 
+ 
+ 
+ 


17, 

time 

time  = 

lx,e9 

lx,e9 


]  ') 


(proutv ( 
,el2.5,6 
3)/19x,6 
3)/19x,6 
lx,e9.3)/19x,6 
lx,e9.3)/19x,6 
lx,e9.3)/19x,6 
lx,e9.3)/19x,6 
lx,e9.3)/19x,6 
lx,e9.3)/19x,6 
lx,e9.3)/19x,6 
lx,e9.3)/19x,6 
lx,e9.3)/19x,6 
lx,e9.3)/19x,6 
lx,e9.3)/19x,6 
lx,e9.3) ) 


i),i=l, 
(lx,e9. 
(Ix,e9. 
(Ix,e9. 
(Ix,e9. 
(Ix,e9. 
(Ix,e9. 
(Ix,e9. 
(Ix,e9. 
(Ix,e9. 
(Ix,e9. 
(Ix,e9. 
(Ix,e9. 
(Ix,e9. 
(Ix,e9. 


numout ) 
3)/19x, 
3)/19x, 
3) /19x, 
3)/19x, 
3)/19x, 
3)/19x, 
3)/19x, 
3)/19x, 
3)/19x, 
3)/19x, 
3)/19x, 
3)/19x, 
3)/19x, 
3)/19x, 


6(lx, 
6(lx, 
6(lx, 
6(lx, 
6(lx, 
6(lx, 
6(lx, 
6(lx, 
6(lx, 
6(lx, 
6(lx, 
6(lx, 
6(lx, 
6(lx, 


e9.3 
e9.3 
e9.3 
e9.3 
e9.3 
e9.3 
e9.3 
e9.3 
e9.3 
e9.3 
e9.3 
e9.3 
e9.3 
e9.3 


)/ 

)/19x, 

)/19x, 

)/19x, 

)/19x, 

)/19x, 

)/19x, 

)/19x, 

)/19x, 

)/19x, 

)/19x, 

)/19x, 

)/19x, 

)/19x, 


6(lx,e9. 
6(lx,e9, 
6(lx,e9, 
6(lx,e9, 
6(lx,e9. 
6(lx,e9, 
6(lx,e9. 
6(lx,e9. 
6(lx,e9, 
6(lx,e9, 
6(lx,e9, 
6(lx,e9, 
6(lx,e9. 


3)/ 

3)/ 
3)/ 
3)/ 
3)/ 
3)/ 
3)/ 
3)/ 
3)/ 
3)/ 
3)/ 
3)/ 
3)/ 


return 

end 

c 
************************************************************************ 

*  * 

*  read  in  node  &  element  data  ,  material  properties  * 

*  ace.  &  force  data,  initial  condition  data,  and  output  request   * 

*  * 

************************************************************************ 

c 

subroutine  readata  ( nnd,nel,nummat, numout, iacc,ndof, mode, 
+  i force, imesh,xc,yc, zc,rifix,e, rkout, inital) 


implicit  real*8  (a-h,o-z) 

dimension  rifix(l) , mode (12, 1) , rkout (3,1) , xc (1) ,yc (1) 

dimension  ifixx(3) ,node (12) ,kout (3) ,e(ll,  1) 


zc(l) 


common  /box  1/  iprob, delta, alpha, toler, gravity 

common  /box  4/  npnts,numif ,nnaf ,ndisi,nveli, 
i-  kfpnts(6)  ,  jnode(30)  ,  jdir(30)  ,  jaf  (30)  , 

i-  ndnod(10)  ,kdis(10)  ,nvnod(10)  ,kvel(10) 

common  /box  4a/g,disi (10) ,veli (10) , f (10) , omega (10)  , 
h  ff (500,6) ,ta(500,6) ,tf (500,6) ,aa(500,6) 
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meq  =  nnd*ndof 
c 

do  100  i=l,meq 
100  rifix(i)  =  0.0 
c 

c  read  &  write  nodal  coord.  &  B.C. 
c   skipped  nodes  will  be  automatically  generated 
c 

write(6,550) 

1  =  0 
110  read(5,*)   n,xc (n) ,yc (n) , zc (n) , (ifixx(i) , i=l, 3) 
c 

do  120  j=l,3 

nl  =  (n-l)*ndof 

rifix(nl+j)  =  dble(if ixx( j ) ) 
120  continue 
c 

nl  =  1+1 


c 

c 


if  (1  .eq.  0)  go  to  130 

dl  =  n-1 

dxl  =  (xc(n)-xc(l) )/dl 
dyl  =  (yc(n)-yc(l))/dl 
dzl  =  (zc(n)-zc(l) )/dl 
c 

130  continue 

I  =  1+1 

if  (n-1)  180,170,150 
c 

150  xc(l)  =  xc(l-l)+dxl 

yc(l)  =  yc(l-l)+dyl 

zc(l)  =  zc(l-l)+dzl 
c 

do  155  j=l,3 

II  =  (l-l)*ndof 
12  =  (l-2)*ndof 
rifix(ll+j)  =  rifix(12+j) 

155  continue 
c 

if  (imesh  .eq.  1)  then 
write(6,570)  l,xc(l) ,yc(l) ,zc(l) , (int (rifix(i) ) , i=ll+l, 11+3 ) 
endif 
go  to  130 
c 
170   continue 

if  (imesh  .eq.  1)  then 
write (6, 570)  n,xc (n) ,yc (n) , zc (n) , (int (rifix(k) ) ,k=nl+l,nl+3) 
endif 

if  (nnd-n)  180,190,110 
c 

180  continue 

write(6,580)  n 
stop 
c 

190  continue 
c 
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c 

c 


read  &  write  element  connectivity 

write(6,600) 
1  =  0 
200  read (5, *)  n,node(l) , node (2) ,  node (3) ,node (4) , node (5) ,node(6) , 
+         node (7) , node (8) , node (9) , node (10) , node (11) , node (12) ,knt 


do  212  i=l,12 
rnode(i,n)  =  dble (node (i) ) 
212  continue 


245 


250 


255 


nl  =  1+1 
if  (knt 
continue 
1  =  1+1 
if  (n-1) 
mode  ( 1 , 
mode  ( 2 , 
mode  ( 3 , 
mode  ( 4 , 
mode  ( 5 , 
mode  ( 6 , 
mode  ( 7 , 
mode  ( 8 , 
mode  ( 9 , 
mode  ( 1 0 
mode  (11 
mode  (12 
go  to  24 
continue 
if  (imes 
write (6, 


eq.    0)    knt 


260,255,250 

1)  =  mode  (1,1-1)    + 

1)  =   mode  (2, 1-1)    + 

1)  =  rnode(3,l-l)    + 

1)  =   rnode(4,l-l)    + 

1)  =   rnode(5,l-l)    + 

1)  =   mode  (6, 1-1)    + 

1)  =  rnode(7,l-l)    + 

1)  =   mode  (8, 1-1)    + 

1)  =  mode  (9, 1-1) 

,1)  =   rnode(10,l-l) 

,1)  =   mode  (11, 1-1) 

,1)  =   rnode(12,l-l) 
5 


dble 
dble 
dble 
dble 
dble 
dble 
dble 
dble 


(knt) 
(knt) 
(knt) 
(knt) 
(knt) 
(knt) 
(knt) 
(knt) 


+  dble (knt) 


+ 
+ 
+ 


endif 
if  (nel- 


h  .eq.  1)  then 

620)  (k,int(rnode(l,k)  )  ,  int  (mode  (2,  k)  )  ,  int  (mode  (3  ,  k)  )  , 
int  (mode  (4,k)  )  ,  int  (mode  (5,k)  )  ,  int  (mode  (6,k)  )  , 
int  (mode  (7, k)  )  ,  int  (mode  (8, k) )  ,  int  (mode  (9, k)  )  , 
int  (mode  (10,  k) )  ,  int  (mode  (11, k)  )  ,  int  (mode  (12,  k)  ) 
k  =  nl,n) 

n)  260,270,200 


c 
c 
c 


c 
c 
c 
c 


260  continue 

write(6,630)  n 

stop 
270  continue 

read  &  write  material  properties 

do  300  i=l,nummat 

read (5,*)   k,e(l, k) , e (2,k) ,e (3, k) , e (4,k) ,e (5,k) 
read (5,*)   e (6,k) ,e (7, k) , e (8,k)  ,e (9, k)  ,  e (10, k) ,e (11, k) 
write(6,960) 

write (6, 970)  k,int(e(l,k) ) , e (2,k)  ,  e (3,k) , e (4,k) ,e (5, k) 
write(6,980) 

write (6, 990)  e (6,k) , e (7,k) ,int(e(8,k) ) ,e (9, k) , e (10, k) 
300  continue 

read  &  write  no.  of  acceleration  history,  and  data  of  time-acc.  pairs 
in  each  ace.  history 

if  (iacc  .eq.  0)  go  to  760 
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read ( 5 , * )     nacc , npnts 

read ( 5 , * )     g 

write (6, 1170) 

write  (6, 1180)  naccnpnts 
c 

do  721  i=l,nacc 

read ( 5 , * )   ( ta ( j , i ) , aa ( j , i ) , j  =1 , npnts ) 
721  continue 
c 

760  continue 

if  (iforce  .eq.  0)  go  to  762 

if  (iforce  .ne.  1)  go  to  763 
c 

c   read  &  write  no.  of  impact  force  history,  no.  of  nodes,  and  data  of 
c   time-acc.  pairs  in  each  ace.  history 
c 

read (5,*)     numif, nnaf 

write(6,1300) 

write(6, 1310)  numif,nnaf 
c 

if  (numif  .eg.  0)  go  to  768 
c 

do  755  i=l, numif 

read (5,*)   kfpnts(i) 

jf  =  kfpnts(i) 

read (5,*)   (tf ( j , i) , f f ( j , i) , j=l, jf ) 

write(6,1295) 

do  785  j=l,jf 

write(6,1210)  i, j , tf ( j , i) , f f ( j , i) 
785  continue 
755  continue 
c 

768  continue 

write(6,1315) 
c 

c   read  &  write  data  of  position  applied  by  arbitrary  shape  impact  force 
c   function  (  node,  d.o.f.,  history  no.  ) 
c 

do  769  i=l,nnaf 

read (5, *)     jnode (i) , jdir (i) , jaf (i) 

write (6, 1225) jnode (i) , jdir (i) , jaf (i) 
7  69  continue 

go  to  762 
c 

c   read  &  write  data  of  position  applied  by  sinusoidal  impact  force 
c   function  (  node,  d.o.f.,  history  no.  ) 
c 

763  continue 

read (5,*)   nnaf 
c 

do  764  i=l,nnaf 

read ( 5 , * )   j  node (i) , j  dir ( i ) ,f(i) , omega ( i ) 

764  continue 
c 

7  62  continue 

if  (inital  .eq.  0)  go  to  765 

if  (inital  .eq.  2)  go  to  782 
c 

c  read  &  write  data  of  displacement  type  I.e.  (  node,  dof,  I.e. value  ) 
c 
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read (5, *)   ndisi 

write(6,784) 

do  786  i=l, ndisi 

read (5,*)     ndnod(i) ,kdis (i) ,disi (i) 

write (6,783)  ndnod ( i ) , kdis ( i ) , disi ( i ) 
7  86  continue 
c 

if  (inital  .ne.  3)  go  to  765 
c 

c  read  &  write  data  of  velocity  type  I.e.  (  node,  dof,  I.e. value  ) 
c 

782  read (5,*)   nveli 
write(6,787) 

do  788  i=l, nveli 

read (5,*)     nvnod(i) ,kvel (i) , veli (i) 
write (6, 783)  nvnod(i) , kvel (i) , veli (i) 
788  continue 
c 

c   read  &  write  data  of  output  record  requested  (node,  d-v-a-type,  dof) 
c 

7  65  continue 

if  (numout  .eq.  0)  go  to  780 
c 

write(6,1228) 
write(7,1228) 
c 

do  770  i=l, numout 
read(5,*)  (kout ( j ) , j=l, 3) 
write (6, 1229)  i, (int(kout (j ) ) , j=l,3) 
write (7, 122  9)  i, (int (kout (j) ) , j=l,3) 
do  771  j=l,3 
rkout(j,i)  =  dble (kout ( j ) ) 
771  continue 
770  continue 
7  80  continue 
c 

550  format  (/2x, 'card  4',7x, 'nodal  point  data',/, 

+        7x, 'node  no . ' , 3x, 'x-ord' , 4x, 'y-ord' , 4x, ' z-ord' 
+        6x, ■ ifx' ,5x, ' ify ,6x, 'ifz' ) 
570  format  (10x, i4, 2x, f 6 .2 , 2x, f 6 .2, 5x, f 6 .2 , 3i8) 
580  format  (17x, 'nodal  point  error  n  =',i5) 
600  format (/2x, 'card  5     element  data' /, lx, ' ele.no. ' ,3x, 'Nl ' , 

+         3x, 'N2 ' ,3x, 'N3 ' ,3x, 'N4' ,3x, 'N5' , 3x, 'N6 ' , 3x, ' N7 ' , 3x, ' N8 ' , 
+        2x, 'mat' ,2x, 'row' ,2x, 'col ' , lx, 'E-con. ' ) 
620  format  (lx,  i4, 4x, i4, lx, i4, lx, i4, lx, i4, lx, i4, lx, i4,lx, 

+        i4,lx,i4,2x,i2,lx,i4,lx,i4,4x,il) 
630  format  (17x, 'element  number  error  n  =',i5) 
c 

783  format  (i8, 6x, i6, lOx, f 12 . 5) 

784  format  (5x, 'node ', 5x, 'disl-dir' , 8x, ' initial  disp',/) 
787  format  (5x, 'node ', 5x, 'vel-dir' , 8x, ' initial  vel',/) 

c 

960  format  (/2x, 'card   6  &  7     material  property  data',/,10x, 

+ 'material  material  mass      Youngs    Poisson   tensile' /10x, 
+' group  no.  type  no.   density   modulus    ratio    strength') 

970  format  (llx, i3, 7x, i3 , 3x,el0 .4, lx,el0 .4,2x, f 5 . 3 , 3x,el0 . 4) 

980  format  (19x, ' cohesion   phi     yield     tangent  hardening' , /3 Ox 
+' angle    criterion  modulus     rule') 

990  format  (17x,el0.4, 2x, f 6.2, 5x, i2, 5x,e9 . 4, lx, f 5.3) 
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1170  format  (/,2x, 'card  12 ', 6x, 'prescribed  acceleration  card'/) 

1180  format  (lOx, 'total  no.  of  acceleration  histories  =  ', 

+    i5, /, 10x, ' total  no.  of  time-acc.  pairs  in  each  ace.  history  =', 

+   i5) 

1200  format  (/,2x, 'card  13 ', 6x, 'acceleration-history  card',/, 

+  17x, "pair  no. ' , 8x, 'time' , 9x, 'ace' /) 

1210  format  (20x, 16 , llx, i6, 4x,el2 .4, 4x,el2 . 4) 

1220  format  (/,2x, 'card  14 ', 6x, 'nodal  accele.  information  card',/, 

+  17x, '  node','   dir','   ace'/) 

1225  format  (15x, i5, 8x, i5, 14x, i5) 

1228  format  (/,2x, 'card   21    output  information  card',/,14x, 

+  'seq.     node*      d- (0) , v- (1) ,a- (2) ,    x(l) ,y (2) , z (3) ' , /,38x, 
+  'stress- (3) ' , 6x, 'xy (4) ,yz (5) ,xz (6) ' ) 

1229  format  (12x, i4, 6x, i4, 13x, i4, 16x, i4) 

1295  format  (/2x, 'card  12  &  13     impact  force  history  card',/,17x, 

+  ' force  history  no. ' 5x, 'pair  no. ' , 5x, ' time ' , 9x, ' iforce ' ) 

1300  format  (/2x, 'card  11     prescribed  impact  force') 

1310  format  (20x, 'total  no.  of  impact  force  history  =',i5,/, 

+  20x, 'total  no.  of  nodes  applied  by  impact  force  =',i5) 

1315  format  (/2x, 'card  14     nodal  impact  force  information',/, 

+  15x, 'node  no.    x- (1) ,y- (2) , z- (3)    force  history  no .' ) 
c 

return 
end 
c 
c 
************************************************************************ 

*  * 

*  compute  the  rigid  body  rotation  by  cross  product  between  * 

*  undeformed  and  deformed  vectors  * 
* 

*  * 
************************************************************************ 

c 

subroutine  thetafind  (xcdl,ycdl , zcdl,xcd2 ,ycd2 , zcd2 , 
+  xu,yu, zu,xd,yd, zd, thetaX, thetaY, thetaZ) 

c 

implicit  real*8  (a-h,o-z) 

data  phi/3. 1415926535898d0/ 
c 

xau  =  xu  -  xcdl 

yau  =  yu  -  ycdl 

zau  =  zu  -  zcdl 
c 

xad  =  xd  -  xcd2 

yad  =  yd  -  ycd2 

zad  =  zd  -  zcd2 
c 

c  compute  the  length  of  undeformed  shape 
c 

rLuxy  =  dsqrt (xau*xau  +  yau*yau) 

rLuyz  =  dsqrt (yau*yau  +  zau*zau) 

rLuxz  =  dsqrt (xau*xau  +  zau*zau) 
c 

c  compute  the  length  of  deformed  shape 
c 

rLdxy  =  dsqrt (xad*xad  +  yad*yad) 

rLdyz  =  dsqrt (yad*yad  +  zad*zad) 

rLdxz  =  dsqrt (xad*xad  +  zad*zad) 
c 
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c  compute  the  unit  components  of  undeformed  vectors 

c 

c  xy  plane  =  1 

euxl  =  xau/rLuxy 

euyl  =  yau/rLuxy 
c  yz  plane  =  2 

euy2  =  yau/rLuyz 

euz2  =  zau/rLuyz 
c  xz  plane  =  3 

eux3  =  xau/rLuxz 

euz3  =  zau/rLuxz 
c 

c  compute  the  unit  components  of  deformed  vectors 
c 
c  xy  plane  =  1 

edxl  =  xad/rLdxy 

edyl  =  yad/rLdxy 
c  yz  plane  =  2 

edy2  =  yad/rLdyz 

edz2  =  zad/rLdyz 
c  xz  plane  =  3 

edx3  =  xad/rLdxz 

edz3  =  zad/rLdxz 
c 

c  compute  the  rigid  body  rotation 
c 
c  xy  plane  rotation (thetaZ) 

Cz  =  euxl* edyl  -  edxl* euyl 

thetaZ  =  dasin(Cz) 

if  ( (xu*xd  .It.  O.dO) .and. (yu*yd  .It.  O.dO))  then 
if  (thetaZ  .ge.  O.dO)  thetaZ  =   phi  -  thetaZ 
if  (thetaZ  .It.  O.dO)  thetaZ  =  -phi  -  thetaZ 

endif 
c  yz  plane  rotation (thetaX) 

Cx  =  euy2*edz2  -  edy2*euz2 

thetaX  =  dasin(Cx) 

if  ( (yu*yd  .It.  0 .d0) .and. (zu*zd  .It.  O.dO))  then 
if  (thetaX  .ge.  O.dO)  thetaX  =  phi  -  thetaX 
if  (thetaX  .It.  O.dO)  thetaX  =  -phi  -  thetaX 

endif 
c  xz  plane  rotation (thetaX) 

Cy  =  euy2*edz2  -  edy2*euz2 

thetaY  =  dasin(Cy) 

if  ( (xu*xd  .It.  O.OdO) .and. (zu*zd  .It.  O.OdO))  then 
if  (thetaY  .ge.  O.dO)  thetaX  =  phi  -  thetaY 
if  (thetaY  .It.  O.dO)  thetaX  =  -phi  -  thetaY 

endif 
c 

return 

end 
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Appendix  2 

Verification  of  the  three  dimensional  finite  element  code 


Problem  1.    A  rectangular  plate  (or  half  space)  of  elastic  material  subjected  to  ramp 
loadings.  Deflections  are  compared  with  the  solutions  obtained  by  using  ANSYS. 

Problem  2.     A  rectangular  plate  (or  half  space)  of  elastic-plastic  material  subjected  to 
ramp  loadings  (Mises  criterion  and  associated  flow  rule).  Progress  of  the  plastic  zone  are 
shown.    Deflections  are  compared  with  the  solutions  obtained  by  using  ANSYS. 

Problem  3.     A  rectangular  plate  (or  half  space)  of  elastic-plastic  material  subjected  to 
ramp  loadings  (Drucker-Prager  criterion  and  associated  flow  rule).  Progress  of  the 
plastic  zone  are  shown.    Deflections  are  compared  with  the  solutions  obtained  by  using 
ANSYS. 

Problem  4.    A  rectangular  plate  of  elastic-plastic  material  with  Mises  criterion  subjected 
to  sinusoidal  loadings  (Response  is  in  the  elastic  range). 

Problem  5.    A  rectangular  plate  of  elastic-plastic  material  with  Mises  criterion  subjected 
to  pulse  loadings  (Response  is  in  the  elastic  range). 

Problem  6.     A  rectangular  plate  of  viscoelastic  material  of  Maxwell  type  subjected  to 
ramp  loadings 


